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Abstract 


Electrochemical energy storage is needed for many mobile technical systems, 
such as communication and electromobility. Lithium-ion batteries (LIBs) have 
attracted intensive attention for electrochemical energy storage over the past 
decades. The enormous utilization of LIBs combined with the limited and un- 
evenly distributed lithium source will drive up the price of lithium. In contrast 
to lithium, sodium-ion batteries (NIBs), which are based on the wide availabil- 
ity, abundance and low cost of sodium, become a potential promising alterna- 
tive to LIBs. In most storage materials, phase changes occur during regular 
operation and, thus, cannot be avoided. The respective phases of a storage 
material possess different lattice constants giving rise to a strain mismatch of 
up to a few hundred percent which, in turn, causes mechanical stresses and, 
thus, leads to damage of the electrode particles. In view of this, it is the objec- 
tive of this work to study phase changes and mechanical stresses in electrode 
particles by means of phase-field simulations. Thermodynamical phase-field 
models rely on continuous order parameters, thus, leading to diffuse interfaces 
between adjacent phases with no need for the cumbersome tracking of the po- 
sition of a sharp interface. 

In this work, for non-linear Cahn-Hilliard type models describing diffusion, 
a thermodynamical framework for the coupling with mechanics is presented. 
For the mechanical part, besides small deformation theory two different finite 
deformation elasticity formulations based on elastic Green strain and logarith- 
mic elastic strain, respectively, are introduced and compared. First, a phase- 
field model for the cathode material Na,FePO, of NIBs is studied for the first 


Abstract 


time to understand phase changes, mechanical deformation, and stress evolu- 
tion in Na,FePOy electrode particles. As a major novelty, the key parameters 
in the phase-field model for Na, FePO. are determined. In this way, our model 
captures the important feature that distinguishes Na,FePO, from Li,FePO,. 
Furthermore, we study the particle size and average concentration dependent 
miscibility gap of the nanoscale insertion materials LiyMn204, Li,FePOa, and 
Na,FePO4. Finally, we introduce the nonlocal species concentration theory 
in terms of diffusion and phase changes from a nonlocal free energy density, 
compare this theory to the Cahn-Hilliard theory, and show how the nonlocality 


influences the results. 
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Kurzfassung 


Für viele mobile technische Systeme wie Kommunikation und Elektromobilität 
werden elektrochemische Energiespeicher benótigt. Lithium-Ionen-Batterien 
(LIBs) haben in den letzten Jahrzehnten intensive Aufmerksamkeit für die elek- 
trochemische Energiespeicherung auf sich gezogen. Die enorme Nachfrage an 
LIBs in Kombination mit der begrenzten und ungleichmäßig verteilten Verfüg- 
barkeit von Lithium wird den Lithiumpreis in die Hóhe treiben. Im Gegen- 
satz zu LIBs werden Natriumionenbatterien (NIBs), die auf der breiten und 
reichlichen Verfügbarkeit sowie den niedrigen Kosten von Natrium beruhen, 
zu einer vielversprechenden Alternative. In den meisten Speichermaterialien 
treten Phasenänderungen während des regulären Betriebs auf und können nicht 
vermieden werden. Die jeweiligen Phasen eines Speichermaterials besitzen 
unterschiedliche Gitterkonstanten, die zu einer Dehnungsfehleranpassung von 
bis zu einigen hundert Prozent führen. Dadurch entstehen mechanische Span- 
nungen, die zu einer Schádigung der Elektrodenteilchen führt. Vor diesem 
Hintergrund ist es das Ziel dieser Arbeit, Phasenänderungen und mechanische 
Spannungen in Elektrodenteilchen mittels Phasenfeldsimulationen zu unter- 
suchen. Thermodynamische Phasenfeldmodelle basieren auf kontinuierlichen 
Ordnungsparametern, was zu diffusen Grenzflächen zwischen benachbarten 
Phasen führt, ohne dass die Position mit einer scharfen Grenzfläche aufwändig 
verfolgt werden muss. 

In dieser Arbeit wird für nichtlineare Cahn-Hilliard-Modelle, die die Diffu- 
sion beschreiben, ein thermodynamisches Gerüst für die Kopplung mit der 
Mechanik vorgestellt. Für den mechanischen Teil werden neben der Elastiz- 


itatstheorie der kleinen Verzerrungen auch zwei verschiedene Formulierungen 
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der Elastizität für große Verzerrungen, basierend auf dem Green'schen sowie 
dem logarithmischen Verzerrungsmaß, verglichen. Zum ersten Mal wird ein 
Phasenfeldmodell für das Kathodenmaterial Na,FePO4 von NIBs untersucht 
um Phasenänderungen, mechanische Deformation und Spannungsentwicklung 
in Na,FePO, -Elektrodenteilchen zu verstehen. Als große Neuheit werden die 
Schlüsselparameter im Phasenfeldmodell für Na,FePO4 bestimmt. Auf diese 
Weise erfasst dieses Modell die wichtige Unterscheidung zwischen Na, FePO4 
und Li,FePO4. Darüber hinaus wird die Partikelgröße und die durchschnit- 
tliche konzentrationsabhángige Mischungslücke der nanoskaligen Insertions- 
materialien Li,Mn> O4, Li,FePO4 und Na,FePO, untersucht. Abschließend 
wird die nichtlokale Spezieskonzentrationstheorie in Bezug auf Diffusion und 
Phasenänderungen aus einer nichtlokalen freien Energiedichte eingeführt, mit 
der Cahn-Hilliard-Theorie verglichen und gezeigt wie die Nichtlokalität die 


Ergebnisse beeinflusst. 
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1 Introduction 


1.1 Generalities on battery systems 


Energy conversion and storage have become a growing global concern over the 
past decades due to urgent energy demand. The renewable energy sources, e.g., 
solar radiation, wind, and waves, are becoming increasingly popular, although 
they can be available in intermittent manner [1]. To mitigate this adversity, 
the chemical energy may help. The chemical energy is the most appropriate 
form of energy storage in terms of energy density, and batteries provide stored 
chemical energy with the ability to deliver it as electrical energy with high con- 
version efficiency and no gaseous exhaust [2]. 

Depending on the rechargeability, batteries are classified into two categories: 
the primary and secondary batteries. Primary batteries evolve an irreversible 
process, converting the chemical energy into electrical energy only once. On 
the contrary, secondary batteries, also known as rechargeable batteries, are pos- 
sibly used several times for repeated discharges and charges. This is due to the 
reversible electrochemical reactions in which the batteries can also be oper- 
ated in the opposite direction to charge them. Therefore the electrical energy 
is converted back to chemical energy during charging so that the batteries may 
be discharged again. We will discuss the reversible process in Section 1.2. 
However, it should be noticed that even for secondary batteries, capacity fade 
of batteries occurs after the repeated charging and discharging cycles owing 
to unavoidable irreversible processes like side reactions, bulk and interfacial 
resistances of both electrodes and electrolyte, and inhomogenous electrode de- 


formations due to the speices insertion and extraction [3]. Hence, increasing 
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the number of battery cycles is nowadays crucial in order to improve the life 
time of batteries. 

Lithium-ion batteries (LIBs) are widely used and common battery type. As 
illustrated in Fig. 1.1, LIBs exhibit superb performances in terms of energy 
density and design flexibility among the existing battery technologies. Indeed, 
lithium is the lightest metallic element (0.53gcm ^?) and has a very low redox 
potential (E Qu u^ — 3.04V v.s. standard hydrogen electrode), which enables 
batteries with high voltage and high energy density [1]. 
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Figure 1.1: Comparison of various battery technologies in terms of volumetric and gravimetric 
energy density [4]. Image reused with the permission of Springer Nature. 


Since the commercial LIBs were first announced by Sony in 1991, LIBs have 
contributed considerably to applications that require portability, such as enter- 


tainment, computing and telecommunication equipment. Furthermore, driven 


1.2 Design and working principle of lithium-ion batteries 


by the improved desire for green technologies, the use of batteries has ex- 
panded from portable electronics to large scale applications, for example, elec- 
tric vehicles [5], which would reduce the pollution and secure energy indepen- 


dence. 


1.2 Design and working principle of 
lithium-ion batteries 


A battery is composed of several electrochemical cells, connected in series 
and/or in parallel, to achieve the required voltage and capacity [6]. As depicted 
in Fig. 1.2, a LIB cell consists of current collectors, anode, cathode, and sep- 
arator between them. A separator that is commonly made of polymers, in an 
ideal case, is an ionic conductor and electronic insulator. The ionic conduction 
of the separator is provided by its pores which are filled with the liquid elec- 
trolyte [7]. A current collector, made of copper or aluminum, contacts with the 


electrode and conducts the electrons to the external circuit. 


5 «——6& — 
discharge charge 


Carbon black 
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positive current collector 
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Electrolyte 
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p 4 Active particle 
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Porous Separator 

Figure 1.2: (Left) Schematic of a LIB cell, consisting of current collectors, cathode, anode, and 
separator. (Right) Illustration of the lithium reduction at the surface of the active par- 
ticles. The lithium-ions are transported through the electrolyte, which fills the pores in 
the separator and electrodes. Electronic conduction occurs in the conductive network, 
made of carbon black and active particles [8]. 
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The electrodes are usually made of a composite of different materials, as shown 
in Fig. 1.2. The reacting material is formed as small active particles. These 
mono- or polycrystalline active particles effectively store the lithium-ions [9, 
10], and are embedded in the soft and porous binders. A binder, usually made 
of polyvinylidene difluoride (PVDF), allows to transfer electrolyte through its 
pores and provides mechanical integrity [11]. To maintain a high electronic 
conductivity of the electrode, the binder is enriched by so called carbon black, 
which is a carbon additive of nanosize [7]. An electrically conductive network 
is built by active particles and carbon additive through which the active par- 
ticles are connected to the current collector. Insertion materials are usually 
chosen as the active material. The structure of this material does not change 
during lithium insertion. As a result, this material acts as a host to provide 
empty lattice sites that can be occupied by the lithium. For the cathode ma- 
terials, like lithium manganese oxide Li;Mn5O4 (LMO), lithium cobalt oxide 
Li,CoO», lithium nickel oxide Li; NiO», and lithium iron phosphate Li,FePO4 
(LiFPO), all of them possess their own merits and demerits [9, 10, 12-15]. It 
should be noticed that LiFPO does not belong to the category of metal ox- 
ides. LiFPO has a lower electrochemical potential than the metal oxides, but is 
thermally more stable [10]. Carbon-based anodes are still very popular anode 
material in the market today. Alloy anodes such as silicon, aluminum, anti- 
mony, and tin are very promising candidates for the next generation of LIBs 
due to their high theoretical energy density, safe operation potentials, and rela- 
tively low cost [7, 16]. 

The fundamental working principle of a LIB cell is the discrepancy in the 
electrochemical potential of lithium in two electrodes [17]. Lithium diffuses 
from the electrode with a high electrochemical potential of lithium to the elec- 
trode with a low electrochemical potential of lithium. When the battery is 
discharged, the difference in the electrochemical potential of lithium of the 
two electrodes drives lithium-ion to diffuse out of the anode, through the elec- 
trolyte, and into the cathode. To keep the electrodes electrically neutral, elec- 


trons flow through an external circuit from the anode to the cathode. Both the 
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ionic and the electronic processes are reversed when the battery is charged by 


an external power source. 


1.3 Motivation and goals 


The interplay between electrochemistry and mechanics in the battery system 
has attracted growing scientific interest in the recent years. During charge or 
discharge, the active material suffers from inhomogeneous volume change, and 
as a consequence, stresses are induced in the active particles of the electrode. 
The stresses may lead to fracture of electrode particles. For example, in LIB 
systems a volume expansion of more than 300% of silicon particles has been 
observed [18], which can cause particle fracture. Examples for the observation 
of fracture in storage particles for different materials are shown in Fig. 1.3. 
Such mechanical degradation leads to capacity loss of a battery over charge 
and discharge cycles [3]. On the other hand, for thermodynamical reasons, 
there is a contribution of the stresses to the driving force for diffusion in the 
host material [19, 20]. 

Intercalation electrode materials commonly exhibit phase changes during in- 
sertion of the intercalation species, such as crystalline silicon (Si) [21], tin, 
antimony and their oxides [18, 22] for anodes, as well as cathode materials 
LMO [23] and LiFPO [24-26]. Phase changes may induce large concentration 
gradients at a mesoscopic scale and thus also large stress magnitudes. Even at 
low C-rates, large stresses may also be induced in the phase-separating elec- 
trode materials [27, 28]. For the modeling of phase changes, it is very common 
to use the Cahn-Hilliard theory where the order parameter is a conserved quan- 
tity in order to avoid the need for the cumbersome tracking of the position of a 
sharp interface [29]. We will review the derivation of the Cahn-Hilliard theory 


from different approaches in chapter 4. 
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Figure 1.3: (a) SiSn thin film is cracked after the first few cycles of lithiation and delithiation 
[30]. Image reused with the permission of Electrochemical Society, Inc. (b) A LiFPO 
particle is cracked into parts after 60 cycles [31]. Image reused with the permission of 
Elsevier. 


LIBs have attracted intensive attention for electrochemical energy storage over 
the past decades. The massive use of LIBs combined with the limited and 
unevenly distributed lithium source will dramatically increase the prices of 
lithium, and high cost remains a critical problem for the development of LIBs 
[2]. In contrast to lithium, the wide availability, abundance and low cost of 
sodium on earth make sodium-ion batteries (NIBs) suitable for large scale en- 
ergy storage devices in which high energy density becomes less critical [32]. 
For comparison, properties of lithium and sodium are summarized in Table 
1.1. Recently, NIBs have been considered as a promising alternative to LIBs, 
since sodium and lithium exhibit similar chemical properties so that sodium 
chemistry could be applied to a similar battery system, and the fundamental 
principles of the NIBs and LIBs are identical [32]. Just like LIBs, the electro- 
chemical processes in the electrodes of NIBs are also coupled to mechanical 
properties. 
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Category Lithium Sodium Ref. 

Crustal abun- | 20 23,600 Ref. [33] 

dance (ppm) 

Distribution 70% South | Everywhere Ref. [34] 
America 

Cost ($/ton), car- | 5000 150 Ref. [35] 

bonates 

Cation radius | 76 106 Ref. [35] 

(pm) 

Atomic weight | 6.9 23 Ref. [35] 

(g/mol) 

Coordination Octahedral and | Octahedral and | Ref. [35] 

preference tetrahedral prismatic 


Table 1.1: Sodium versus Lithium characteristics. 


Among all cathode materials for NIBs, phosphate based cathode materials are 
among the best candidates, due to their thermal stability and higher voltage 
[2]. Olivine Na,FePO4 (NaFPO) has the highest theoretical specific capacity 
(154mAhg~') compared to the other phosphate based materials (NaV PO4F, 
Na3V2(PO4)2F3 and Na5FePO,4F etc.) [2], which makes it a promising cath- 
ode material for NIBs. Similar to olivine LiFPO, olivine NaFPO also shows 
phase changes during sodium insertion or extraction [32, 36-39]. 

As far as we know, no work has been published for the phase-field modeling 
of the cathode materials of NIBs by now. In this work, a phase-field model for 
NaFPO of NIBs will be studied for the first time. As a major novelty, the mate- 
rial parameters for NaFPO, in particular, the key parameters in the phase-field 
model, are determined. In this way, our model captures the important fea- 
ture that distinguishes NaFPO from LiFPO. Furthermore, the volume expan- 
sion of FePO4 upon sodiation to NaFePO, reaches about 17%, which is much 
larger than that for LiFPO (about 6.8%) changing from FePO4 to LiFePO4 
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[37]. Thus, it is also interesting to study whether or not the small deformation 
theory (SDT) still has a sufficient capacity to represent the deformation of this 
material. 

On the other hand, the reduction of the electrode particle size down to the 
nanoscale range can improve the power density [40, 41] and the rate capability 
[42, 43]. It is found that the miscibility gap between the species-poor phase 
and the species-rich phase shrinks as the particle size decreases, correspond- 
ing to an increase in the mutual solid solubility, and, as a consequence, the 
misfit strain is reduced [44-51]. The thermodynamics of phase segregation 
in nanoscale particles is distinctly different from bulk materials, and the par- 
ticle size dependent miscibility gap plays a nontrivial role in the performance 
of such nanoscale insertion materials. The previous research works about the 
particle size miscibility gap have been almost exclusively focused on LiFPO. 
By now, the particle size and average concentration dependent miscibility gap 
of the nanoscale insertion materials LMO and NaFPO are still unclear. The 
question is still open how the mechanical stress affects the particle size and 
average concentration dependent miscibility gap in intercalation electrode ma- 
terials. In this work, we will investigate and compare the particle size and 
average concentration dependent miscibility gap for the three cathode mate- 
rials LMO, LiFPO, and NaPPO during insertion, using a coupled phase-field 
model based on the Cahn-Hilliard theory and finite deformation elasticity. We 
will physically explain the average concentration dependent miscibility gap, 
and determine the critical particle size below which phase segregation is inhib- 
ited for the three cathode materials. We will also study how the mechanical 
stress affects the particle size and average concentration dependent miscibility 
gap. 

According to Cahn and Hilliard [29], the Cahn-Hilliard theory is derived by 
truncating the Taylor series of the system free energy, ignoring fourth-order 
and higher derivatives of the order parameter in the system free energy. There- 
fore, the system free energy in the Cahn-Hilliard theory depends on the value 


of the order parameter at a certain position and its immediate neighborhood. 
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As a result, the Cahn-Hilliard theory can not be seen as arising from an exact 
macroscopic description of microscopic models of interacting particles [52]. 
In addition, the Cahn-Hilliard equation involves a fourth-order spatial deriva- 
tive of the order parameter, leading to computational difficulty. In this work, 
we will introduce a phase-field theory in terms of diffusion and phase changes 
from a nonlocal free energy density, which can be applied to electrode materi- 
als. This theory is named nonlocal species concentration theory (NSC theory). 
As a special feature, we discuss the nonlocal nature of NSC theory, i.e. in 
which respect NSC theory accounts for effects in the whole problem domain. 
Furthermore, the results from this theory and the Cahn-Hilliard theory are com- 


pared. 


1.4 Overview 


The goal of this thesis 1s to understand phase changes, mechanical deforma- 
tion, and stress evolution in NaFPO electrode particles of NIBs, and to study 
the particle size and average concentration dependent miscibility gap for the 
three cathode materials LMO, LiFPO, and NaPPO. Also, we introduce the NSC 
theory from a nonlocal free energy density, and show how the nonlocality in- 
fluences the results. 

The work is organized as follows. A literature review about related concepts, 
methods, and previous studies is given in chapter 2. In chapter 3, fundamentals 
of continuum mechanics are briefly reviewed. In chapter 4, the derivation of the 
Cahn-Hilliard theory from various approaches is summarized. Furthermore, 
the NSC theory is introduced from a nonlocal free energy density. Finally, we 
discuss the relationship between the CH theory and our theory. In chapter 5, the 
coupled phase-field model of the Cahn-Hilliard theory and finite deformation 
elasticity is derived. For the mechanical part, besides SDT two different finite 
deformation elasticity formulations based on elastic Green strain and logarith- 


mic elastic strain, respectively, are introduced and compared. In chapter 6, a 
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phase-field model for the cathode material NaFPO of NIBs is studied for the 
first time, and the relevant material parameters of NaFPO are determined. We 
present and discuss phase changes and stresses in spherical NaFPO electrode 
particles, and compare the cathode material NaFPO of NIBs to LiFPO of LIBs 
in terms of phase changes and stresses. In chapter 7, we investigate the parti- 
cle size and average concentration dependent miscibility gap of three cathode 
materials LMO, LiFPO, and NaPPO. In chapter 8, we present the results for 
spherical LMO electrode particles based on the NSC theory, discuss the inter- 
face evolution and the nonlocal effect, and compare the results from the NSC 


theory and the Cahn-Hilliard theory. We conclude our work in chapter 9. 
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2.1 Phase-field method 


The phase-field method is a powerful computational approach for modeling 
microstructure evolution in materials at the mesoscale without explicitly track- 
ing the positions of interfaces. It was initiated in 1894 by van der Waals to 
model a liquid-gas system using a density function that varies continuously 
at the liquid-gas interface [53]. Approximately sixty years ago, Ginzburg and 
Landau [54] formulated a phase-field model for superconductivity using a com- 
plex valued order parameter and its gradients, and later Cahn and Hilliard [29] 
carried out a thermodynamic formulation considering the gradients in thermo- 
dynamic properties in heterogeneous systems with diffuse interfaces. 

Microstructures of physical systems are usually compositional and structural 
heterogeneous by nature, and can evolve with time. The driving force for mi- 
crostructural evolution is the possibility to reduce the free energy of the sys- 
tem. The microstructural evolution is described by means of the field variables 
called order parameters that are continuous functions of time and spatial co- 
ordinates. The field variables have different values in different phases to dis- 
tinguish different microstructures. There are two types of field variables: con- 
served and nonconserved. The former one has to satisfy the local conservation 
condition [55, 56]. The temporal and spatial evolutions of the conserved and 
nonconserved field variables are governed by the well-known Cahn-Hilliard 
nonlinear diffusion equation as in the present case and the Allen-Cahn relax- 
ation equation, respectively. Both equations can be physically derived from a 


so-called local microforce balance of Gurtin [57]. The field variables can be 
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an existing physical order parameters such as concentration [29], polarization 
[58], and magnetization [59]. On the other hand, for the sole purpose of avoid- 
ing tracking the interfaces and identifying different phases, dummy fields are 
introduced as field variables. For example, in the cases of solidification [60] 
or crack propagation [61], the phase-field that equals zero represents the liquid 
phase or the unbroken phase while the phase-field with value of unity repre- 
sents the solid phase or the crack phase. 

Phase-field model is based on a diffuse-interface description. The interfaces 
between phases are continuous across the interfacial regions, as shown in Fig. 
2.1a, which is different from the sharp-interface model of the conventional 
approach in Fig. 2.1b. As a result, an impressive advantage of the phase-field 
method is that there is no need to follow explicitly the position of the interfaces 


through mathematical equations during microstructural evolution. 
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Figure 2.1: (a) Diffuse interface: continuously between phases. (b) Sharp interface: discontinuous 
at the interface [55]. Image reused with the permission of Elsevier. 


The application of the phase-field method has been extensively developed in 
the past decades, besides solidification [62, 63] and solid-state phase trans- 
formations [56, 64], the phase-field method also has been employed in grain 
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growth and coarsening [65-67], the spinodal decomposition of binary mix- 
ture [29, 68], crack propagation [61, 69—74], electromigration [75, 76], vesicle 
membranes in biological applications [77, 78], and multicomponent interdiffu- 
sion [79]. 


2.2 Overview of phase-field modeling of 
electrode materials 


Although multi-physical models have extensively investigated in the past decades, 
as shown in comprehensive overviews [80-84], it remains a big challenge to 
simulate the complicated microstructure of electrode materials. The interac- 
tions among mechanical, electrochemical, electrical, and thermal fields in the 
battery system induce complex non linear partial differential equations for the 
unknown variables (i.e. displacement, concentration, electric potential, and 
temperature). As a result, it needs more advanced numerical methods to solve 
such complex multi-physical models. As shown in Figs. 2.2 and 2.3, differ- 
ent processes that occur in a battery during normal operation are described by 
their corresponding mathematical models. In this work, we will focus on diffu- 
sion, mechanical deformation, and phase segregation in active particles at the 
mesoscale, corresponding to processes of number 3, 4, and 5 in Figs. 2.2 and 
2.3, respectively. In what follows, the overview of the phase-field modeling of 


electrode particles will be presented. 
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Figure 2.2: Different processes that occur in a battery during normal operation [84]. Image reused 
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with the permission of Springer Nature. 
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Figure 2.3: Corresponding models for the processes that occur in a battery during normal operation 
[84]. Image reused with the permission of Springer Nature. 
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2.2.1 Phase-field modeling of phase changes 


By now, most of the research works about phase-field modeling of phase 
changes have been focused on LiFPO of LIBs. In 2004, Han et al. [85] used 
the classical Cahn-Hilliard phase-field model to describe phase segregation in 
LiFPO electrode particles, and investigated to what extent non-Fickian behav- 
ior can affect results from experimental techniques for measuring diffusion 
coefficients, such as Galvanostatic Intermittent Titration Technique (GITT) 
and Potentiostatic Intermittent Titration Technique (PITT). According to their 
results, GITT and PITT, based on the Fick's equation, are still capable of accu- 
rately determining diffusion coefficients even if the solid is characterized by a 
gradient energy term. 

In 2008, Singh et al. [86] developed a general continuum theory for phase- 
transformation dynamics in LiFPO electrode particles by coupling an 
anisotropic Cahn-Hilliard phase-field model with surface reactions governing 
the flux of ions across the electrode/electrolyte interface. The general model 
can induce different transport and phase separation dynamics regimes. For 
the regime of isotropic bulk-transport-limited phase transformation dynamics, 
bulk diffusion in all directions is much slower than surface reactions, and the 
phase boundary is fully contained within the material and moves along the 
direction of the ionic flux. As shown in Fig. 2.4 a, the shrinking core structure 
forms. For the regime of anisotropic bulk-transport-limited phase transforma- 
tion dynamics, as depicted in Fig. 2.4 b, diffusion in the x and z directions is 
negligible, while diffusion in the more accessible y direction is much slower 
than surface reactions. The phase boundary still shows a shrinking core in the 
bulk, but ions are confined to 1D channels in the y direction. Thus, anisotropy 
alters the shape of the phase interface rather than its basic diffusive dynamics. 
For the new regime of anisotropic surface-reaction-limited dynamics, surface 
reactions are much slower than diffusion in the y direction but much faster 
than diffusion in the x and z directions. As illustrated in Fig. 2.4 c, the phase 


boundary extends from surface to surface along planes of fast ionic diffusion, 
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rather than forming a classical shrinking core. It should be noticed that these 
three reaction-diffusion arguments are oversimplified and only give a sense of 
possible dynamical regimes in the general model. 


(a) ? 


1 


Figure 2.4: Transport models obtained in different limits of the characteristic timescales for bulk 
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diffusion and surface reactions. Figures show xy cross sections of spherical and plate- 
like single crystals during Li insertion, after phase nucleation has occurred. Lithiated 
portions of the crystal are shaded, and points outside particles represent flux of Li 
ions across the electrode/electrolyte interface (shown only for spherical particles). The 
FePOy/LiFePO, phase boundary is denoted by the dashed line, and arrows indicate 
movement of the boundary as insertion proceeds. (a) Isotropic bulk transport limited. 
(b) Anisotropic bulk transport limited. (c) Anisotropic surface reaction limited [86]. 
Image reused with the permission of Elsevier. 
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Subsequently, Bai et al. [87] developed a Cahn-Hilliard phase-field model of 
reaction-limited intercalation in anisotropic LiFPO nanoparticles and predicted 
that phase separation is suppressed above a critical current. As shown in Fig. 
2.5, spinodal decomposition or nucleation leads to moving phase boundaries 
for small currents. Above a critical current density, the spinodal decomposi- 
tion is found to disappear, and the particles start to fill homogeneously, which 


may explain the superior rate capability and long cycle life of nano-LiFPO 
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Figure 2.5: Voltage responses at different constant currents [87]. Image reused with the permission 
of American Chemical Society. 


Recently, Santoki et al. [88] used the classical Cahn-Hilliard equation to study 
the mesoscopic effect of the surface curvature of the cathodic particle made of 
LMO. They found that, near the convex region of the particle surface, more 


sites are available for the applied flux than for the sites that host lithium ions. 


17 


2 Literature review 


As a result, the onset of phase segregation prefers to occur in high curvature 
regions of the particle. Also, the elliptical particle with a higher aspect ratio is 


subjected to the onset of the phase segregation, prior to the lower ones. 


2.2.2 Phase-field modeling of diffusion and 
mechanics 


Here, a related overview of phase-field modeling of diffusion and mechanics 
is given (see also Zhang and Kamlah [89]). A general phase-field theory for 
the coupling of the Cahn-Hilliard equation and finite deformation elasticity 
based on the local dissipation inequality and the so-called local microforce 
balance has been presented by Gurtin [57]. In this theory, the external micro- 
force itself remains a quantity free to choose, and by taking it equal to zero, the 
classical Cahn-Hilliard equation is obtained. Walk et al. [28] have identified 
the influence of the external microforce in the chemical potential for spherical 
LMO particles. Cogswell and Bazant [50] investigated the effects of elastic 
coherency strain on the thermodynamics, kinetics, and morphology of interca- 
lation in LiFPO nanoparticles based on a reaction-limited phase-field model. 
The key parameters in the phase-field model for LiFPO are determined. Their 
model quantitatively captures the influence of misfit strain on solubility seen 
in experiments. They concluded that the elastic coherency strain strongly sup- 
presses phase separation during discharge, enhancing the rate capability and 
extending the cycle life. Tang et al. [90] used a Cahn-Hilliard phase-field 
model wih SDT to investigate phase separation and crystalline-to-amorphous 
transformations in spherical isotropic LiFPO nanoparticles, and assessed the 
conditions under which amorphous phase transitions may occur in LiFPO par- 
ticles. Huttin and Kamlah [27] considered the coupling between the Cahn- 
Hilliard equation and SDT for spherical particles of LMO, and demonstrated 


that large stresses may also occur even at low C-rates. 
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Anand [91] derived a theory for lithium insertion modeled by the Cahn-Hilliard 
equation combined with finite deformation elasto-plasticity based on the multi- 
plicative decomposition of the deformation gradient, in which the external mi- 
croforce is identified as related to the Mandel stress in the chemical potential. 
To this end, the so-called macroforce and microforce balances are evaluated by 
following the virtual-power approach of Germain [92] and Gurtin [93]. Subse- 
quently, Di Leo et al. [94] formulated a continuum model which coupled the 
Cahn-Hilliard-type phase-field theory with finite deformation elasticity based 
on logarithmic elastic strain. However, Anand's model [91, 94] ignored the 
fact that logarithmic elastic strain is dependent on the species concentration in 
the stress chemical potential, although logarithmic elastic strain is also a con- 
stitutive variable. As a result, the extremely complicated and important term 
induced by the concentration dependence of the logarithmic elastic strain in 
the stress chemical potential is not accounted for. The concentration depen- 
dence of the strain tensor in the chemical potential is also confirmed by Gurtin 
[57]. Interestingly, the external microforce term related to the Mandel stress 
in Anand's model [91, 94] is just equivalent to the ignored term induced by 
the concentration dependence of the logarithmic elastic strain, which also is 
verified by Refs. [28] and [95]. Walk et al. [28] reported that for SDT the 
identification of the external microforce in Anand's model [91, 94] leads to a 
doubling of the mechanical coupling term in the chemical potential. Zhao et al. 
[95] derived a Cahn-Hilliard phase-field model coupled to mechanics based on 
the Neo-Hookean elasticity without introduction of the external microforce in 
the chemical potential, and showed that their model agrees with Anand's model 
[91, 94] which ignored the concentration dependence of the strain tensor in the 
chemical potential. As remarked by Gurtin [57], no constitutive relation is 
specified for the external force. Rather, specifying the external microforce by 
invoking the so-called principle of virtual power may be considered an option. 
However, in view of the above discussion we do not find it appropriate to adopt 


this in our work. 
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2.2.3 Studies of the miscibility gap in nanoscale 
insertion materials 


Here, a related overview of the studies of the miscibility gap in nanoscale in- 
sertion materials is given (see also Zhang and Kamlah [96]). The miscibility 
gap is the difference between the lower and upper bounds of the concentration 
range of phase segregated states. The research works about the particle size 
dependent miscibility gap have been almost exclusively focused on LiFPO of 
LIBs. Meethong et al. [45] observed that the miscibility gap of LiFPO shrinks 
as the particle size decreases, and they suggested an estimate of the critical par- 
ticle size (15 nm) below which a complete solid solution is achievable at room 
temperature for LiFPO. Wagemaker et al. [49] revealed that the miscibility gap 
for LiFPO in small particles not only shrinks, but also strongly depends on the 
average concentration, which is also called "state of charge" (SOC). In con- 
trast to our common thermodynamic knowledge that the solubility limits are 
independent of the average concentration, their combined neutron and X-ray 
diffraction investigation reveals an decreasing miscibility gap that appears to be 
strongly dependent on the average concentration below particle sizes of 35 nm. 
However, they explained the average concentration dependent solubility limits 
based on the so-called average solubility limit (i.e. the average compositions in 
each phase). In our understanding, the average concentration dependent mis- 
cibility gap should be related to the minimization of system free energy. What 
is more, the influence of the elastic strain energy is ignored in their diffuse in- 
terface model. The elastic strain energy can suppress phase segregation, and 
cause larger solid-solution-composition-ranges [27, 50]. Cogswell and Bazant 
[50] used a coupled phase-field model based on SDT to study the particle size 
dependent solubility for LiFPO. They performed the phase-field simulations by 
allowing a square particle at SOC- 0.5 to relax to equilibrium at zero current. 
Not looking at other SOC values, the average concentration dependent misci- 


bility gap is not accounted for. Welland et al. [51] developed a comprehensive 
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phase-field model based on SDT accounting for facet dependent surface wet- 
ting and investigated equilibrium states of LiFPO particles of different size and 
average lithium concentration, and found that the miscibility gap vanishes for 
particles of a radius around 5 nm, and that the solubility limits change with 
overall particle lithiation in small particles. However, no flux is applied at the 
surface of the particles in their work, rather a phase segregated concentration 
profile or a fluctuation field was used as an initial condition for different sur- 
face wetting cases. They have not yet addressed the dynamic loading, which 
is related to the experimentally more relevant condition of a constant applied 
flux. Zhao et al. [97] developed a phase-field model accounting for lithium 
transport in particles, phase separation, and interface reactions across the parti- 
cle network to study the influence of particle size variation on the performance 
of LIBs with V205 nanowires. Their work reveals that both particle size and 
size variation in electrodes should be small to avoid intra- and inter-particle 


phase separation. 


2.2.4 Nonlocality of phase-field models 


Here, a related overview of nonlocality of phase-field models is given (see also 
Zhang and Kamlah [98]). Starting with Walter [99], nonlocal models have 
been studied in the papers [100-111]. In the above mentioned works [99— 
107, 109-111], a basic form of the nonlocal energy is used to describe the 
behavior of materials that exhibit a morphology of phase structures. Compared 
to the second-order derivative of the order parameter in the system free energy 
in the Cahn-Hilliard theory, the nonlocal energy form in the nonlocal model 
has no derivative of the order parameter so that both, discontinuites and oscil- 
lations, are allowed. The nonlocal energy is written in an integral form. To 
be specific, it is related to a weighted integral average of the squared devia- 
tion between the order parameter at the point of consideration and the order 


parameter anywhere else in the problem domain. This will be introduced in 
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Section 4.2.1. In this way, as the system free energy contains an additional 
term representing an energy average over the entire problem domain, the con- 
cept of nonlocality is introduced. The nonlocality is based on the consideration 
of nonlocal (long-range) particle-particle interactions due to molecular forces 
from the viewpoint of both classical and modern statistical mechanics [99]. 
Fosdick and Mason [112] derived a continuum stored energy that includes a 
nonlocal consideration based on the work of Brandon [103]. 

Urbachs et al. [113] derived a nonlocal diffuse interface model for the mi- 
crostructure evolution of tin-lead solder, in which a nonlocal mass fraction is 
used in order to introduce the nonlocal effect. However, the integral equation 
related to the normalized weight function in their work [113] does not satisfy 
a natural weight function property [103, 112, 114]. We will discuss the nat- 
ural weight function property in Section 4.2.1. As a result, the definition of 
the nonlocal mass fraction seems to be not physically sound. Indeed, by their 
expression of what they call the nonlocal mass fraction the nonlocal species 
concentration, in our wording is not equal to the species concentration for ho- 
mogenous states, if their theory satisfies the natural weight function property. 
Also, the nonlocal mass fraction should be explained as a weighted average 
value of the local mass fraction in the entire problem domain rather than, as 
they say, over a finite zone in the vicinity of the interface. What is more, 
the definition of the interface tension coefficient seems to be mathematically 
inconsistent with the other definitions. Di Leo et al. [94] formulated the gra- 
dient micromorphic concentration theory using the principle of virtual power, 
in which besides conventional species concentration an additional micromor- 
phic concentration is employed. However, this theory is not derived from the 
nonlocal approach mentioned above, and the relationship between their theory 
and the Cahn-Hilliard theory seems to not to be clear. Therefore, they just 
pursued the goal how the solutions from their theory can approach the solu- 
tions from the Cahn-Hilliard theory. In particular, they do not discuss aspects 
of the nonlocal nature of their theory which can account for effects beyond an 


infinitesimal neighborhood of the point of consideration. 
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mechanics 


3.1 Kinematics 


For the description of deformation, a body ZZ is considered in the Euclidean 
space. A reference configuration is arbitrarily chosen at a fixed reference time. 
Usually, an undeformed body is preferred to be the reference configuration, 
and this also corresponds to the initial configuration at time t = fo. X denotes 
an arbitrary material point of Z in the reference configuration described by 
Lagrangian or material coordinates, and X is the spatial position of this material 
point in the current configuration described by Euler or spatial coordinates. A 
smooth one-to-one mapping from the material points X to the spatial point s X 


at time f describes the motion of Z 
X — X(X,t). (3.1) 
The displacement field 
i —x—X (3.2) 


represents the connection vector between the current and reference configura- 


tions. The displacement gradient H is expressed by 


H = Gradi. (3.3) 


23 


3 Basics of continuum mechanics 


The differentiation of equation (3.1) leads to the deformation gradient 
F = Grad Y = H + I, (3.4) 


which is very crucial in nonlinear continuum mechanics. Here Grad is the 
gradient operator calculated with respect to the reference configuration, and I is 
the second order unit tensor. The deformation gradient F is a two-point tensor 
involving points in two different configurations, and it is a primary measure of 
deformation [115]. As showin in Fig. 3.1, the fundamental relation between 


the material line element dX and the spatial line element dX 
dx =F dX (3.5) 


is determined by the deformation gradient. 


Reference configuration Current configuration 


Figure 3.1: Transformation of a material line element from the reference configuration to the cur- 
rent configuration. 


The determinant of the deformation gradient J gives the mapping from the 


infinitesimal reference volume dVr to the current state dV 


dV = JdVg = detF dVr. (3.6) 
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By means of equation (3.6) the following relation is obtained 
dV — dA. dX = JdÁg- dX, (3.7) 


with dA = dA ii and dÁR = dAr rig denoting the infinitesimal surfaces in the 
current and reference configurations, respectively. Here, n and nr are the cor- 
responding the unit normal of the surfaces. 

Using Equations (3.5) and (3.7) the mapping from the infinitesimal reference 


area dÁR to the current area dA can be obtained as 


dA = JF ! dA. (3.8) 


3.2 Stress tensors 
For infinitesimal surfaces, 
df =p dA = P dAg (3.9) 


holds for the infinitesimal force df, where the tractions p and P are defined in 


the current and reference configurations, respectively, read as 
p=Tn and P— Tg rig. (3.10) 


Here, Tg is the second-order first Piola-Kirchhoff stress tensor, which char- 
acterizes the force per unit undeformed area acting on surfaces in the current 
configuration, and T denotes the symmetric Cauchy or true stress tensor. 
Combining Equations (3.9) and (3.10) yields 


T dÄ — Tg dÁg. (3.11) 
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Using Equation (3.8), the relation between the two stress tensors can be rewrit- 


ten as 
TR=JTF'. (3.12) 


As shown in Equation (3.12), Tr is not symmetric. Using a pullback operator 


F-!, the symmetric second Piola-Kirchhoff stress tensor T¢ is introduced 
T =F T. = JF TF, (3.13) 


and T° describes the force per unit undeformed area acting on surfaces in the 
reference configuration. 

In small deformation framework, the stretch is small, the rotation is close to 
unity, and the difference between material and spatial are neglected, i.e the 


displacement gradient H satisfies 
[H| « 1, (3.14) 


and |-| denotes the norm of a vector. The linear strain tensor is defined as the 


symmetric part of the displacement gradient 


£= 1+0). (3.15) 
3.3 Mechanical equilibrium 


An arbitrary body 2 with the boundary 0.2 is subjected to the body force 
b, and the traction vector is distributed over its surface. The body is in me- 


chanical equilibrium when the resultant force is zero. The following equations 


26 


3.3 Mechanical equilibrium 


describe the equilibrium conditions in the spatial and material presentations, 


respectively, 

[vias | bav =o. (3.16) 

OB B 
n Ty iin dAg | have 0. (3.17) 

BR BR 

Using the Gaussian integral theorem follows 
[ (av +b) av =0, (3.18) 

B 


J, Div Te +5) av, =0. (3.19) 
4 R 


Since the considered body is arbitrary, the above relations are always supposed 
to hold. Hence, the specified local forms of the above mechanical equilibrium 


conditions can be deduced 


dvT+b=0, (3.20) 


Div Ta +b — 0. (3.21) 
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4 Phase-field theory of phase 
changes 


4.1 Cahn-Hilliard theory 


The Cahn-Hilliard theory can be physically derived from the thermodynamics 
by two main approaches. One is based on a variational formulation, follow- 
ing Cahn and Hilliard [29, 68]. The other is derived from the so-called local 


microforce balance [57]. Both approaches will be reviewed in the following. 


4.1.1 System free energy 


To start with, the system free energy is described (see also Zhang and Kamlah 
[89]). We introduce as an order parameter, depending continuously on space, 
the species concentration c, which is measured in mol per unit volume. The 


free energy density consisting of two parts is given by 


> 


y (E, Vc) = v"? (c) + wf" (Ve), (4.1) 


where c is the dimensionless concentration, which is normalized with respect 
to the maximum species concentration Cmax as € = C/Cmax, and V denotes the 
Nabla operator. The homogeneous free energy density y"? is a multiwell 


potential defining the respective phases. Based on the references [27, 85], the 
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homogeneous free energy density exhibits a double-well structure. The exis- 
tence of this zone of concavity indicates that homogeneous species concentra- 
tion states do not always ensure the system free energy to be minimal. In the 
concavity zone, the system becomes unstable towards phase segregation. The 
Maxwell construction, which connects the neighborhoods of the two minima 
by a common tangent, predicts that the system splits into a two phase system 
and the chemical potentials of both the species-poor phase and species-rich 
phase are equal. In this sense, the homogeneous free energy density is a mul- 
tiwell potential where the wells define the respective phases. The second term 
on the right hand side of Equation (4.1) represents the gradient energy leading 


to a diffuse interface between two adjacent phases, which is given by 
oio l, 2 
yë (Ve) = kpTrefNaCmax(54|Vel J: (4.2) 


Here kg is the Boltzmann constant, N4 is the Avogadro constant, and T,.y is the 
reference temperature. Furthermore, À is the gradient energy coefficient with 
units of length squared controlling the thickness of the diffuse interface. 


The system free energy of an arbitrary body & of volume V is 
we) = | qv") wi (We))av. (43) 
B 


The dimensionless free energy density is introduced as V = Y /(kgT,efNACmax), 
and any parts of y may be normalized in the same way. According to Cahn and 
Hilliard [29], the system free energy (c) can be interpreted physically such 
that, to the first approximation, the free energy of a small volume of nonuni- 
form solution can be expressed as the sum of two terms, one being the free 
energy that this volume would have in a homogeneous solution and the other a 


“gradient energy". 
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The homogeneous free energy density for a two phase material can be ex- 
pressed as [27, 85] 


[07 
w""P(c) = kpTserNACmax (a+ De 
T 
+— (ēInē+ (1—c)In (1 2). (4.4) 
Tref 


The first two terms on the right hand side of Equation (4.4) represent the inter- 
action energy, where positive values of o characterize the energy of inserting 
a species into the host material, and negative values of œ indicate the inter- 
action of neighboring species to be attractive. Concerning 1, it is hard to 
relate a physical meaning to non-positive values of oj. On the other hand, it 
has to be noted that a; does not occur in the Cahn-Hilliard evolution equation 
for concentration. According to Huttin [116], there is a critical temperature, 
namely T, = —1/402T,ef, above which the homogeneous free energy density 
is of convex shape, representing an ideal solution. For temperatures T below 
the critical temperature, the free energy density exhibits a zone of concavity 
where the homogeneous concentration states are not stable states of the sys- 
tem, and phase segregation becomes possible, as shown in Fig. 4.1. The zone 
of concavity corresponds to the condition that the second order derivative of 
the homogeneous free energy density with respect to the concentration is neg- 
ative. This inequality is never fulfilled if œ is not negative. Thus, for a system 
of noninteracting species, i.e. OQ = 0, or for a system of species that repel each 
other, i.e. Q > 0, the homogeneous free energy density keeps a convex shape 
for all T > 0. For an attractively interacting species system, i.e. Q2 < 0, de- 
pending on the temperature, the homogeneous free energy density may exhibit 
a zone of concavity. Therefore, at T = T,.y, the attraction is strong enough 
to initiate phase segregation for oo < —4. The terms multiplied by absolute 
temperature T represent the entropy of mixing, where cInc represents the en- 


tropy of mixing valid at low concentration, and (1— c)In(1 — c) describes the 
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non-dilute solution behavior and represents the entropy of mixing responsible 
for the saturation effect [95, 117]. 


—— T2150 K 
—— T=T =298.15 K 
—— T=T,=387.595 K 
—— T2450 K 


imp 


Figure 4.1: Normalized homogeneous free energy density versus normalized concentration for 
different values of absolute temperature T. 


4.1.2 Derivation from a variational formulation 


The standard derivation of Cahn-Hilliard theory starts from the system free 
energy, as defined in Equation (4.3), and the variation Ó with respect to the 


order parameter c reads 


Q yp kpT.>¢N, m 5 

GP = " Me dg bv a) dV 
B dc Cmax 

| ow kpTretNa 

B 


Oc Cmax 


AV?c) óc dV 


keTrefNa „> 
«f Bref A Yen Sc dA. (4.5) 
OB 


Cmax 
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Here V? is the Laplacian operator. Equation (4.5) yields the natural boundary 


condition on the entire boundary 0.4 
Ve-n=0. (4.6) 


When the interface between phases meets the particle surface, Equation (4.6) 
enforces that it is perpendicular to the surface [95]. Thus, the chemical poten- 
tial is defined as the variational derivative of the total free energy with respect 
to the species concentration 
hit ð mwp k T, N 
u=% = ká B ref A y Y2, (4.7) 


c Qc Cmax 


and it represents a thermodynamic driving force for diffusion and phase 


changes. Thus, Equation (4.5) can be rewritten in this form 


ay™? — kpTofN. 
Jy Y B ref A) 26) 8c dV 
B 


dc Cmax 


ow 


25 u óc dV. (4.8) 


Based on the conservation of mass, the evolution of species transport is given 
by 


é = — div(J) (4.9) 


with the mass flux J related to the chemical potential u through the Onsager 


relation 


J-2-M-Vy, (4.10) 
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where the mobility tensor M is non-negative definite. we choose an isotropic 


mobility according to 
M(c)=M (c)I (4.11) 


with the function 


= Doc (Cmax = c) 


M(c)= 
( ) kpT,efNACmax ' 


(4.12) 
which is symmetric in the range between zero and maximum concentration. 
Here, Do is the diffusion constant. 

Combining Equations (4.7), (4.9), (4.10), and (4.11) yields the classical Cahn- 


Hilliard equation 


x mwp 
é — div (mov (= keTresNa avec) ) (4.13) 


oc Cmax 


4.1.3 Derivation from a microforce balance 


According to Gurtin [57], although the above derivation of Cahn-Hilliard equa- 
tion is simple and physically sound, it should not be regarded as basic, but 
rather as precursor of more complete theories. Indeed, such derivation limits 
the manner in which rate terms enter the equations, and requires a-priori spec- 
ification of the constitutive equation (4.3). Here, we will review the work of 
Gurtin [57] to show that how to derive the Cahn-Hilliard equation from a mi- 


croforce balance. 
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In the Cahn-Hilliard theory the kinematical process is related to the order pa- 
rameter c, and it is reasonable that there should be “microforces” whose work- 
ing (expenditure of power) accompanies changes in c. Given an arbitrary con- 
trol volume R (subregion of a body 2), the system of forces is presumed to be 


consistent with the microforce balance [57] 
farna] É ii dA —0, (4.14) 
R OR 


where & is a vector stress, and 7 and y represent, respectively, internal and 
external forces distributed over the volume of a body 2. 
The above global law (4.14) should be satisfied for all time and all control 


volumes R, and the local microforce balance is obtained 
divé - x4 y -0. (4.15) 


The general Cahn-Hilliard equation can be derived based on balance of mass, 
the local microforce balance, and a generalization of the dissipation inequality 
accommodating diffusion. In order to be consistent with the balance of mass 
(4.9), here we ignore the external mass supply. Actually, the external mass 
supply term will be automatically dropped out in the following local dissipation 
inequality even if it is taken into account. 

According to the second law of thermodynamics, the rate at which the free 
energy of R increases can not be more than the working on R plus the rate at 


which free energy is carried into R by mass transport. The working is given by 


wq) - [ve av e | &- aa. (4.16) 


and the chemical potential through the contributions 


- [ufa (4.17) 
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characterizes the rate at which free energy is carried into R by mass transport. 
Therefore, the appropriate form of the second law of thermodynamics is the 


dissipation inequality 


[wavs [rave f € 8)edA-- the MEN (4.18) 


The above dissipation inequality can be manipulated further 


fiw- frea- fE aeda | LJ fü dA 
= [wav f veav - | div (eB) av + | div (uJ) dV 
R R R R 
& ] vv - frav- [ (e iv £- 8-2) av 
R R R 
+ [in div J4- J- Vu) dV 
R 


< 0. (4.19) 


Using the local microforce balance (4.15) and the balance of mass (4.9), Equa- 
tion (4.19) yields the local dissipation inequality 


y -é-Ve-J-Vu 4 (x— u)é x 0. (4.20) 


According to Gurtin [57], the constitutive equations are considered of the form 


= E(c,Ve,u, Vu), x — f(c, Vc, u, Vp). (4.21) 


All the above constitutive equations should be consistent with the second law 
of thermodynamics in the form of the local dissipation inequality (4.20) , which 


leads to the following inequality 


E J-Vp <0, (4.22) 
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holding for all fields c and u. As a result, we can induce the following relations, 


which are sufficient for the above inequality (4.22), 


> QU(c, Vc 
ft(c, Vc, u) = u m vo , 
EN See i (c, V 
E = êle, Ve) = VOM) (4.23) 


where y and ë are independent of u andVu, f is independent of Vu, and 
J.Vuxo (4.24) 


for all fields, such that J' has the form of Equation (4.10). Using Equations 
(4.23) and (4.15) yields 


dlc, Vc) 
Ve 


|— Oc, Vc) 
B oc 


u div( )- Y. (4.25) 
It can be interestingly seen that Equation (4.25) with y — 0 gives the chemical 
potential obtained from the variation derivative of the total free energy with 
respect to species concentration, as shown in Equation (4.7). 

Substituting Equations (4.10) and(4.25) into the balance of mass (4.9) yields 


the generalized Cahn-Hilliard equation 


é — div (mov (sut GE, JJ (4.26) 
OVc 


oc 


and the standard Cahn-Hilliard equation is obtained if y = 0. 
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4.2 Nonlocal species concentration theory 


Here, a nonlocal species concentration theory for diffusion and phase changes 
is introduced from a nonlocal free energy density (see also Zhang and Kam- 
lah [98]). This theory can be interpreted as an extension of the Cahn-Hilliard 
theory. In principle, nonlocal effects beyond an infinitesimal neighborhood are 
taken into account. We first introduce a general formulation of a nonlocal free 


energy density, and then derive the NSC theory from the nonlocal model. 


4.2.1 Nonlocal free energy density 


In the nonlocal model, nonlocal interactions are introduced by the weighted 
spatial average of the local concentration. Based on [52, 103, 112], the free 


energy density including truly nonlocal consideration reads as 


y = y" (e(x)) + JA iR o(¥3X) (EX) — e(y))^av. (4.27) 


Here, y"? is the homogeneous free energy density, which has been expressed 
in Equation (4.4). The second term on the right hand side of Equation (4.27) 
is the nonlocal free energy density. Z (Z C IR?) is the entire problem domain 
containing all material points with position vectors X. In the integral, y denotes 
the spatial variable for the volume integration over the entire problem domain 
B. Furthermore, A is a material coefficient for the nonlocal particle-particle 
interaction penalty with units J /m? [112] and &(Y;X) is a positive, symmetric 
weight function that weights the relative energetic contribution of concentra- 
tion fluctuations to the free energy density. The nonlocal free energy density in 
Equation (4.27) vanishes for homogeneous states, i.e. when the concentration 
at all positions is the same. Thus, the so-called nonlocality becomes active in, 
say, a phase segregated state, i.e. when the concentration at any position y, 


possibly far away from X, is different from the concentration there. The weight 
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function @(¥;x): IR? — R* is a monotonically decreasing function of the dis- 
tance p between X and y. It assumes recognizably non-vanishing values over 
a finite zone in the neighborhood of the material point X, only, and asymp- 
totically approaches zero as the distance p grows to infinity. In this way, the 
weight function attributes less weight as the distance p increases. This satisfies 
the physical idea that the ability of a material point y to influence the local free 
energy at a material point X decreases as the distance p between the two mate- 
rial points increases. Therefore, the weight function possesses a natural weight 


function property [103, 112, 114] 
o(y;X)dV — 1, (4.28) 
R3 


where, again, y is the variable of spatial integration. According to Brandon 


[103], the weight function &(y; X) ensures that the operator 
ec L'(Z;R)o 1 ox(y;x)e(y)dV e LP (Z;R) (4.29) 
is compact, where the space L^(4; R) is the set of functions © 
c: BR (4.30) 


such that 


1/4 
Il € lloc) Ue] < o. (4.31) 


A similar definition can be applied to the space p (428;R). For example, the 
type of the weight function can be chosen of the form [102, 103, 112, 114] 


o;z) ~ (4.32) 
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which satisfies the natural weight function property and compactness condition 
when n > 0. The weight function (4.32) is only dependent on the distance 
p = |y-X]. Taking the variation of the free energy (4.27) with respect to c(X) 


to obtain the chemical potential then gives 


SW gy"? A "ee 
bon = Eom um t1], 903068)-:094v, (433) 
óc(X) Oc(X) Cmar JB 


where, similar as in Equation (4.3), 


we L v(y)aV (4.34) 


is the system free energy of the entire problem domain 2. 


We now define two parameters according to 


- A 
a=, | 06DAN, (4.35) 
CmaxB B 
A zt 
= " o: X)av. (4.36) 
Cmax JB 


Here é = €/Cmax is the normalized nonlocal species concentration that is a 
weighted average value of c in the entire problem domain Z. The term 
A/(cmaxß) scales the normalized nonlocal species concentration © such that 
© is equal to the normalized species concentration © for homogenous states. 
The parameter f is a penalty energy coefficient with units J/mol, which is 
related to the penalty energy density, which will be introduced later. It can be 
normalized as B = p / (kpT; NA). 

Using Equations (4.35) and (4.36), the chemical potential can be rewritten to 
yield 

ov oy"? 


Pon = Sy = 3e + kpTre ¢NaB (E) — EE). (4.37) 
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In analogy with Peerlings et al. [114], two different partial differential equa- 
tions (PDEs) can be derived from the integral formulation (4.35), which will 


be discussed now in detail. 
4.2.2 Helmholtz equation governing nonlocal 
species concentration 


The integral formulation (4.35) can be written in a differential form. First, c(y) 


is expanded into a Taylor series around X 


26) = Me + LHR): G- 36-3) 
ea EGG -)6—:)6—5) 
ve: G—3)8-3) 8-3) 8-3). (438) 


Again, V denotes the Nabla operator with respect to spatial position X. Sub- 
stituting Equation (4.38) into the integral formulation of the nonlocal species 


concentration (4.35), yields 


ida za i 5h osx) Ve()) G — x)dv 

Er a. que , 6:3) V Ves) : (¥—X)(¥—X)dV 

= ssl )VVVe(3): Gi —3) (F — X) (5—x)av 

i a VVe(3) =: (y-3)(9 - X) (y — x) (95 — X)dV 
n (4.39) 
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If an isotropic and homogeneous weight function is considered, as the example 
in Equation (4.32), which is only dependent on the distance p = |y — X], i.e. 
(¥3X) = o(|y — ¥|) = o (p), odd derivatives vanish upon integration leaving 


&(x) = EE) - PV?e(x) + FV EE) + IG VSe(x) 4---- (4.40) 
with 

1 A 

26 xp od 2 

l^ = 2 Crab n o (p)dV, (4.41) 
1 A 

4. inc VE 4 

b = Al css nr o (p)dV, (4.42) 
1 A 

Éb = / $ 4.4 

3 Lo p @(p)dV, (4.43) 


22222 


so on. Since the series in Equation (4.40) is not truncated, its validity is not 
limited to some infinitesimal neighborhood of position X, in principle. In this 
sense, nonlocal species concentration € is truly nonlocal. 

Next, we apply the Laplacian operator to Equation (4.40) and multiply with /? 


to obtain 
PVPE) = PV?c(x) + V EE + PVE EZ) +... (4.44) 
Subtracting Equation (4.44) from Equation (4.40) gives 


E- PVE) = + — I )v^e(x) 
u$ PB) +... (4.45) 
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For vanishing coefficients of the higher order derivatives of c, Equation (4.45) 


reduces to the inhomogenous Helmholtz equation [114] 
ĉ- PV =. (4.46) 


In Appendix A.1, we discuss under which conditions the inhomogenous Helm- 
holtz equation (4.46) becomes the governing equation of nonlocal species con- 
centration ©. To this end, we first show that the Green's function of Equation 
(4.46) satisfies the natural weight function property (4.28). Since additionally 
the Green's function is isotropic and homogeneous, it is chosen as the weight 
function for nonlocal species concentration ©. Second, we show that ĉ defined 
in this way satisfies the inhomogenous Helmholtz equation (4.46), indeed. 
According to Equation (4.41), the parameter / is of the dimension length, and 
in view of its role in Equation (4.46), we call it characteristic interface length 
scale. It is taken to be on the order of magnitude of the interface thickness d 
between adjacent phases of the species concentration field c [113], and mea- 
sures the volume in the neighborhood of the material point X where species 
concentration c contributes significantly to the nonlocal species concentration. 
For example, for a weight function of the type in Equation (4.32), by Equa- 
tion (4.41) there is a relation between / and the parameter 7. Therefore, the 
characteristic interface length scale / is related to the scale of microstructure. 
Henceforth, we assume / to be a constant. The Helmholtz equation (4.46) is a 
PDE to calculate the nonlocal species concentration. The Helmholtz equation 
is an elliptic partial differential which represents a time-independent form of 
the wave equation obtained by separation of variables. It arises, for example, 
in the study of electromagnetic radiation, seismology, and acoustics. 

In order to solve the Helmholtz equation (4.46), we need to introduce a bound- 


ary condition. Here the natural boundary condition 


Vé.i-0 (4.47) 
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is used with ri being the unit normal to 0.2. For an interpretation of this bound- 


ary condition, integration of Equation (4.46) over the domain Z gives 


| édV — i, PV? ĉdV = / cdV. (4.48) 


The second term of Equation (4.48) vanishes by the divergence theorem and 


boundary condition (4.47), which gives 


cdV = / édV. (4.49) 


Thus, the natural boundary condition (4.47) preserves the global species con- 
tent of the entire problem domain in the nonlocal averaging [114]. In other 
words, at each time instant, the total content of nonlocal species concentration 
é is equal to the content of the considered species in the problem domain 2. 

In summary of the above derivation, together with the boundary condition 
(4.47), the partial differential equation (4.46) for ĉ is an alternative to the inte- 
gral representation (4.35) when using the Green's function as weight function. 
Therefore, the NSC theory based on the Helmholtz equation (4.46) is a special 
form of the nonlocal model, in which the weight function o (y;X) is defined as 
O(Y;X) = G(y;X). Indeed, the infinite series of higher-order derivatives of c is 
still implicitly present in the gradient term on the left-hand side of Equation 
(4.46) [114]. As mentioned before, this indicates that spatial interactions in- 
duced by the higher-order derivatives of c represent effects at a finite distance 


and so the NSC theory is truly nonlocal. 


4.2 Nonlocal species concentration theory 


4.2.3 Decomposition of nonlocal free energy 
density 


The system free energy 


y = wen 


= f, y" (c )+ yperaliy (c(3),66)) a yrs (&(y))av 
(4.50) 


of the NSC theory consists of three parts. For a derivation see Appendix A2. 
The first part of Equation (4.50) is the free energy that the entire problem do- 
main would have in a homogenous solution. In the second part, 
yrenalty (2(%),2(&)) is the penalty energy density of the form 
x 1 - 

y?" (s), EE) = 5 emab (ER) EM)". (4.51) 
This part considers a system penalty energy induced by the the difference be- 
tween normalized species concentration c(X) and normalized nonlocal species 


concentration x). A similar penalty term is also introduced in [94] and [118]. 


For the last part, we introduce a variance energy density of the form 


range (E) Tenab 1 o(p —é(X 3))^qv. (4.52) 


The last part of V accounts for a system variance energy which depends on the 
weighted average value of the squared deviation of the local species concentra- 
tion c(y) from normalized nonlocal species concentration &(X). The variance 
energy density is related to the concept of variance in the probability theory 
and statistics. It should be noticed that the variance energy density makes no 
contributuion to diffusion and phase changes due to the independence of the 
variance energy density with respect to c(X). Compared to the free energy den- 


sity (4.27) in the nonlocal model, the nonlocal free energy density is split into 


45 


4 Phase-field theory of phase changes 


two terms in the NSC theory, i.e. the penalty energy density y^*"4"* and the 


variance energy density y’#ance, 


4.3 Comparison of two diffusion theories 


For a better comparison between these two diffusion theories, we first derive 
the Cahn-Hilliard theory from the nonlocal model. If the terms of order four 
and higher of Equation (4.40) are neglected, Equation (4.40) yields an approx- 
imation of Equation (4.35): 


6-64 D V. (4.53) 


Using Equation (4.53), Equation (4.37) gives the chemical potential, as shown 
in Equation (4.7), for the Cahn-Hilliard theory with the gradient energy coef- 
ficient A = BP, which is assumed to be a constant, henceforth. According to 
Equation (4.53), the nonlocal species concentration in a material point depends 
only on the species concentration and its second-order derivative. As a conse- 
quence, the spatial interaction induced by only the second-order derivative of 
species concentration ¢ is limited to the infinitesimal neighborhood of the con- 
sidered material point. Thus, the Cahn-Hilliard theory is weakly nonlocal. 

Compared to the NSC theory, the Cahn-Hilliard theory is a simplified version 
of the NSC theory, since the NSC theory implicitly accounts for an infinite se- 
ries of higher-order derivatives of ¢ in Equation (4.46) while the Cahn-Hilliard 
theory only contains the second-order derivative in Equation (4.53). In addi- 
tion, as will be shown in chapter 8, the interface thickness d is controlled by the 
two parameters B and / in the NSC theory, while it depends solely on param- 


eter A in the Cahn-Hilliard theory. Now we compare the corresponding field 
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equations from the two diffusion theories. Combining Equations (4.9), (4.10), 
and (4.37) yields the diffusion equation for the NSC theory as 


3 ay" INS 
ċ = div (mow ( àc T kpTrefNaB e-8)) 3 (4.54) 


Here, Equation (4.54) is a partial differential equation (PDE) which involves 
second-order spatial derivatives in c, and the nonlocal species concentration is 
calculated by solving the Helmholtz equation (4.46). Thus, the NSC theory 
consists of two coupled second-order PDEs. On the other hand, the Cahn- 
Hilliard equation (4.13) is a fourth-order PDE. Therefore, the Cahn-Hilliard 
theory needs stronger continuity requirements on the species concentration 
than the NSC theory. Instead of a fourth-order equation, the two second-order 


PDEs of the NSC theory are computationally less demanding. 
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5 Cahn-Hilliard theory with 
mechanics 


In the previous chapter 4 we introduced the Cahn-Hilliard theory without me- 
chanics. In this chapter, we will derive the coupled model of the Cahn-Hilliard 
theory and mechanics. We first start with the coupled Cahn-Hilliard theory 
with SDT. Then, the coupled Cahn-Hilliard theory with finite deformation elas- 
ticity is derived, including two finite deformation elasticity formulations based 
on elastic Green strain and logarithmic elastic strain, respectively. Besides, for 
the spherically symmetric boundary value problem, the mathematical formula- 


tions of the coupled Cahn-Hilliard theories are also derived. 


5.1 The coupled Cahn-Hilliard theory with 
small deformation theory 


5.1.1 System free energy 


According to Gurtin [57], in order to take the coupling between diffusion and 
mechanics into account, the coupling energy or elastic strain energy should be 


added into the system free energy 


V(c,gradc,£) = n 
B 


] 0O +y (grade) + v^ (e,c)dV. 6.1 
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Here y""? and y8? have been previously referred as the multiwell potential 
and the gradient energy density, respectively. The coupling energy density y^? 


is assumed to be given by 


Cp omi e ARE . e€ 
y 2 £^:C:€ 
= Gle’-e° tr£^)?). 5.2 
(E°: e° s (re*)?) (52) 
Here C is the elasticity tensor 
2Gv 
C=G (54.551 + 515 jx) + [ody Ou Se (5.3) 


which is taken to be constant and isotropic. Accordingly, v is the Poisson’s 


ratio, and G is the shear modulus 


E 
G = ——— 5.4 
2(1+v)’ en 
where E is the Young’s modulus. 
The elastic strain €^ is given by 
£ —£E£—£, (5.5) 
where 
, 1 
E = 39(c- co)I (5.6) 


is the stress-free strain induced by species insertion or extraction, and € is the 
total strain tensor that has already introduced in Equation (3.15). Here, Q is 


the partial molar volume, and co is the initial species concentration. 
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5.1.2 Linear elasticity 


The Cauchy stress tensor can be derived from the free energy density [57], 


which leads to the law of linear elasticity 


T 


d V(c,gradc, €) 


PI 


dw? (£,c) 


dE 


5.1.3 Field equations 


C: £* 


(5.7) 


The chemical potential considering the coupling between diffusion and me- 


chanical deformation is a superposition of three terms 


ow 


hcc cepi (5.8) 
with 
ay" : awn 
mwp 
s dc A (5 — 
x =0 
T 
= kpTretNa (v + OnE 4- T. (Inc — In (1 2)) ; (5.9) 
ref 
gd gd 
psd = oy iv ow: 
dc d gradc 
=0 
= —kpTyepNad div (gradé) , (5.10) 
: oy? ow? 
cp L 
BUT rct Hy (ec) 
e =0 
= cm (5.11) 
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where Ty = 1/3T;; is the hydrostatic stress. Combining Equations (4.9), (4.10) 
and (5.8) yields the coupled Cahn-Hilliard diffusion equation 


é — div (M(c)V(u""? + p + WP), (5.12) 


where M(c) has been defined in Equation (4.12). 
It should be noticed that, replacing the term u8? shown in Equation (5.12) by 


the term 


enalty 

penalty __ dy? ? 

u Lu en, 
óc 


= kpTrepNaB (c — €), (5.13) 
the coupled NSC diffusion equation with SDT can be obtained. 
Finally, based on the balance of momentum, as shown in Equation (3.20), the 


mechanical equilibrium condition is 
div T — 6, (5.14) 


where the body force is neglected. 

Combined with the constitutive equations introduced above, the field equations 
can be taken as partial differential equations for concentration c and displace- 
ment vector u, which need to be solved for given initial and boundary condi- 


tions. This is a fourth-order nonlinear initial-boundary-value problem. 


5.1.4 Mathematical formulation of the spherically 
symmetric boundary value problem 


In this part, the phase segregation problem of species insertion into or extrac- 
tion from a cathodic particle is mathematically formulated. In order to avoid 
costly three-dimensional simulations, here, a spherical particle of radius Ro un- 
der spherically symmetric boundary conditions is considered, and the spherical 


coordinate system (r, 0, 9) is introduced. It is assumed that not only the particle 
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geometry, but also the unknown variables, such as the species concentration c 
and the displacement vector i, as well as the boundary conditions holding at 
the electrode particle surface are invariant under rotation with respect to both 
the 0 and @ coordinates [116]. As a result, the three-dimensional problem is 
allowed to be replaced by an equivalent one-dimensional problem, as sketched 
in Fig. 5.1. Due to the spherical symmetry, all fields are expressed as a function 
of the time ¢ and the radial coordinate 0 < r < Ro: 


c — c(nt), (5.15) 
ii = u,(r,t)é,. (5.16) 
n zl R=Ry 
Spherical symmetry c -— — — À M 
[o o GA] 
c=c(r,t) 


ii — u(r,t)é, 


(Spherical coordinates) 


Figure 5.1: Spherical particle of the spherically symmetry: the three-dimensional problem reduces 
to a one-dimensional problem. 


In numerical implementations, the scaling of the equations by material param- 
eters may lead to larger numerical errors, especially when the prefactors are 
very small or large. Therefore, the boundary value problem is implemented in 
a normalized form. In this work, the relevant normalized parameters are listed 
in Table 5.1. 
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Parameters Normalization 
Concentration c e¢=— 
Cmax 
Radius r F=% 
0 
Radial displacement u, ü, = Re 
Time £ t= 2 
u 0 
Partial molar volume Q Q = Qemax 
I m 
Flux J, J, = mn 


Young's modulus E 


= E 
kg Tref Na Cmax 


Free energy density y 


y — y 
y= KpTye f NACmax 


Chemical potential u 


"ES 


ee e 
kpgTyo p NA 


Table 5.1: Normalized parameters. 


For the spherically symmetric boundary value problem considered here, the 


strain tensor is diagonal, i.e. entries only exist on the main diagonal. The radial 


and tangential components of it are given by (referring to Equation (3.15)) 


E, 


Er 


(5.17) 


(5.18) 


According to Equation (5.5), the elastic strain €* is expressed as 


&£ = 8-8 
ur — 2 (c — eg) 0 
= 0 u& — 9 (c— eg) 0 
0 ur — 9 (c — co) 
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(5.19) 
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Furthermore, based on the law of linear elasticity (5.7), the components of the 


stress tensor are given by 


v 
Te = 26 e y (e 2) 
1 3v 
5 (1+) ae a). 
(5.20) 
V 
T = 2G e+ 7 (e+e) 
1 3v 
5 (1+) ale a) 
(5.21) 


Based on Equation (5.14), the mechanical equilibrium condition of the spheri- 


cally symmetry is formulated as 


a 42r, — m) — Q0. (5.22) 
or r 


Now we derive the coupled diffusion equation of the spherically symmetric 
boundary value problem. For the chemical potential, u”? has been expressed 
in Equation (5.9). According to Equations (5.10) and (5.11), u and u‘? are 


expressed as 


us! = —kpTepNaA div (grad?) 
ae 209 
= kp Trep Nad 5 Tus (5.23) 
u? = —QTy 
2 2 gu, 4.u 
ed Q? Q r Fr 
G É (eere) 3 or 3 r 
2vOQ Qu, Ur 
Q 2 : 24 
ex (c co) + 3 + 23] (5.24) 
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The gradient of the chemical potential is the driving force for diffusion, and the 


components of it are obtained as: 


gradu"? = u Hill e»t. (5.25) 
gradu? = —KpTreyNaA (5: de x)* (5.26) 
vat = o [rt dno alte) 
ea 
(5.27) 


Substituting the above Equations (5.25)-(5.27) into Equation (5.12), we can 


derive the final dimensionless coupled diffusion equation 


TX 3) 9c | 20% 2:96 
pl Nae Fae POF 


óc 2.9, 4. /lou, a 
iz 2 r r Uy 
rel Da ER 3952 a ) 


2vQ 49€ a uy 1 di, Uy 
oa Qn + ap a( IF DE (5.28) 


We can find that the coupled Cahn-Hilliard diffusion equation involves fourth- 


| 
N 


order spatial derivatives in the concentration and third-order spatial derivatives 
in the radical displacement. Here, if the term related to shear modulus G is not 
taken into account, the above coupled diffusion equation (5.28) will degenerate 


to the pure Cahn-Hilliard diffusion equation. 
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N | Pa Jii R=0 du 


u=0 T.i-0 
— — 
dc >, 
—=0 J +n = const. 
Or 
ye | ES He gade- ni =0 
Ors 


Figure 5.2: Boundary and symmetry conditions for the spherically symmetric boundary value 
problem. 


We now discuss the boundary and initial conditions for the spherically symmet- 
ric boundary value problem (see also Zhang and Kamlah [89]). The boundary 
conditions are sketched in Fig. 5.2. At the surface, the particle is assumed to 


be stress free, 
T(Ro,t)-i=0, (5.29) 


where ñi refers to the outgoing unit vector normal to the particle surface. We 


choose a spatially independent mass flux at the surface as 


5 CemaxRo for ¢(Ro,t) <c 
J(Ro,t) -7 = 260053 (Ro,t) S Cmax . (5.30) 
0 for c(Ro,t) = Cmax 


Here C is the C-rate, and C — n means that the amount of a species of a fully 
charged particle would flow into the particle within 1/n hours. Once the maxi- 
mum concentration Cmax is reached at the surface, the mass flux will be stopped. 
Here neglecting surface wetting [119], the natural boundary condition is ex- 


pressed as 


grade-n = Rot) — 0. (5.31) 
or 
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In addition, at the particle center, the boundary conditions 


u, (0,1) — 0, (5.32) 
Oc 
<< (0.1) — 0. (5.33) 
9? 
5,0 -0 (5.34) 


are to be satisfied. Here, Equation (5.32) guarantees the continuity of the parti- 
cle at the center, and Equations (5.32)-(5.34) are needed to ensure the spherical 
symmetry. What is more, these three boundary conditions ensure the physical 
requirement that the flux at the particle center vanishes [116]. 


Finally, the initial conditions are given by 


u, (7,0) — 0, (5.35) 


c (r,0) — 0. (5.36) 


5.2 The coupled Cahn-Hilliard theory with 
finite deformation elasticity 


5.2.1 Kinematics 


The kinematics developed here is related to the basic fields shown in Table 5.2, 
some of which have been discussed in chapter 3.1 (see also Zhang and Kamlah 
[89, 96]). 
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Descriptions Basic fields 
Motion ¥ = %(X,t) 
Displacement vector i—x-X 
Deformation gradient F = Grad X 
Volume ratio J — detF 
Multiplicative decomposi- | F = F*F* 

tion of F 

Elastic deformation gradient | F*,J* = det F° 
Concentration induced de- | F?,J? = det F? 
formation gradient 

Elastic right stretch tensor U‘ = v Fe Fe 
Elastic left stretch tensor Ve = v Fere" 


Polar decomposition of F* 


F° = R*U^* = VRS 


Elastic Green strain tensor 


E° = 5 (F°'F*—I) 


Spectral decomposition of 
U‘ 


Ut 2 Y? | A Fa fg 


Spectral decomposition of 
ve 


v= I Aa (RFa) e (R°Fa) 


Logarithmic elastic strain 


tensor 


e _YV3 ez P 
log = PM In AG Ta @ ra 


Spatial logarithmic elastic 
strain tensor 


ber = Le. Ag (RP) © (Ra) 


Table 5.2: Basic fields in the kinematics. 


Consider a motion X = Y (X ,t) of a homogeneous isotropic body &. The kine- 


matical basis for the coupling between finite deformation elasticity and diffu- 


sion is the multiplicative decomposition [91] 


F — FF, 


(5.37) 
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as illustrated in Fig. 5.3. Here F* is the elastic deformation gradient that is 
a mapping from the intermediate configuration onto the current configuration, 
and the concentration induced deformation gradient F? caused by species in- 
sertion or extraction is a mapping from the reference configuration onto the 


intermediate configuration. F? is assumed to be purely volumetric in the form 
FS = B'I = 4/14 O(cg — co)L (5.38) 


where cg is the species concentration, measured in mol per unit reference vol- 
ume. 
Using Equations (3.6) and (5.37) yields 


J=J 5, (5.39) 


where J* = detF* and J’ = detF$ = 1 4- O(cg — co). 


Intermediate configuration 


m — oy 


p 


F^ 


Reference configuration 
Current configuration 


Figure 5.3: Multiplicative decomposition of F [89]. 
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The right and left polar decompositions of F° are expressed by 
F° = R°US = VR‘, (5.40) 


where R° is a rotation, and U* and V° are the symmetric, positive-definite 


elastic right and left stretch tensors: 
U’ = V Fe! F? and V^ = VFeFe'. (5.41) 
In addition, the right elastic Cauchy-Green tensor is 
C=U?=F'F. (5.42) 


The elastic Green strain tensor is then defined as 


E = (cn 
MP 1 eT me 
- dev) 
= 5 (+ 2(ce—co)) 3 E F1) (5.43) 


Next, the spectral representations of U* and V* are 


3 3 

U = Y Afra Fa and V° = VAs (R74) & (Ro), (5.44) 
a=1 a=1 

respectively where (71,72,73) and (Af, 45, A5) are the orthonormal eigenvectors 

and the positive eigenvalues of U? respectively. According to Anand [91], the 

logarithmic elastic strain tensor and the spatial logarithmic elastic strain tensor 


can be obtained as 


3 
Ej, = InU* =) InAgra & Ta, (5.45) 


a=1 
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e _ e 
Elog.H = InV 


3 
ina; (R75) & (RFq). (5.46) 


a=1 


Using Equations (5.37), (5.38), (5.41), and (5.45), the above logarithmic elastic 


strain tensor can be rewritten as 


= ]nU* 
= InvF'r 


1 
= To) F'F) 


= InVF'F— -In(1-- O(cg — co))I. (5.47) 


e 
log 


Thus, the logarithmic elastic strain tensor involves an eigenstrain accounting 


for the strain developed by composition changes from the reference state. 


5.2.2 System free energy 


The system free energy of the coupled theory, in the reference volume, includes 
three parts [57, 91] 


Vn f Vg d VR 
BR 


A (V? "P (cr) + wi (Grader)+ W^") dV, — (5.48) 
BR 


where yy"? and ysd have been previously discussed as the multiwell potential 
and the gradient energy density, respectively. Here, entities with subscript R 
indicate these to be defined in the reference configuration. The additional term 
V." is the coupling energy density defining the coupling between diffusion and 
mechanics. As shown in Fig. 5.3, the elastic Green strain tensor is defined in 


the intermediate configuration. Thus, the coupling energy of the entire system 
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in terms of elastic Green strain tensor in the intermediate configuration is given 
by 


; 1 
wer = f - E: C: E'dV;. (5.49) 
Bg 2 
According to the multiplicative decomposition of F, we can obtain 
dV, = det F? dVr = J? dVg. (5.50) 


Using Equations (5.49) and (5.50) yields the coupling energy in terms of elastic 


Green strain tensor in the reference configuration 
; 1 
wer zi - JS E° : C: E°dVp. (5.51) 
Ba 2 


Therefore, the coupling energy density in terms of elastic Green strain tensor 


in the reference configuration is given by 


C, 1 S 
wP (cr, E°) = 5 JE:C:E' 
V 
1—2v 


GJ, \E* : Ec + (trEP) |, (5.52) 


Similarly, the coupling energy density in terms of logarithmic elastic strain 


tensor in the reference configuration is derived as 


; 1 
v (cn, Er,,) == 2 J m DC: E; 


" V 2 
= Gs Ej. : Ef. + Tov (trEj,,) 1 (5.53) 
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5.2.3 Elasticity based on elastic Green strain 


According to Gurtin [57], the thermodynamical potential relation for the first 


Piola-Kirchhoff stress Tg is derived from the local dissipation inequality as 


_ OWR 


Te 
A JF 


(5.54) 


The second Piola-Kirchhoff stress T° related to intermediate configuration is 


symmetric considering the symmetry of the Cauchy stress T: 
qT^—J E TRO, (5.55) 
Using Equations (3.12), (5.54) and (5.55), we can obtain the relation 


qe _ 2 9WR _ 1 OWR 
JS aCe JS JEL 


(5.56) 


Using Equations (5.48) and (5.52), Equation (5.56) then yields the isotropic 


elasticity law based on elastic Green strain as 


£ 1 {ower oys (Grader) Pye (cr, E°) 


T° = T 
r| 3E Oe JE? 
—— 
=0 =0 
MELLE 
(Js o Ee 
_ e| M e 
— 2GE 20, (trE^)I 
2 
= 2GE°+ (x- 56) (trE°)I, (5.57) 
where the bulk modulus K is expressed by 
E 
K = —  ——. 5.58 
3(1— 2v) O28) 
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5.2.4 Elasticity based on logarithmic elastic strain 


In the following, we will introduce another isotropic elasticity law based on 
logarithmic elastic strain (see also Zhang and Kamlah [89]). According to 
Anand [91], the second Piola-Kirchhoff stress T° is derived from the free en- 


ergy density according to 


T= 2 OWR 


= ge (5.59) 


Here, Equation (5.59) agrees with Equation (5.56) which is derived from the 
local dissipation inequality of Gurtin [57]. The Mandel stress related to the 


intermediate configuration is defined by 

M = CT“. (5.60) 
Using Equations (5.42) and (5.45) yields 

fog = In(C*?). (5.61) 


Hence, using Equations (5.59) and (5.61), Equation (5.60) then gives 


20WR 1 dyr 
M’ = C° = ; .62 
OF IC T FE, e 


65 


5 Cahn-Hilliard theory with mechanics 


Using Equations (5.48) and (5.53), Equation (5.62) then yields the isotropic 


elasticity law based on logarithmic elastic strain as 


1 |9w7"^ dyk (Grader) IWR (ce Ej,,) 
J | OES, JE: 


log 


1 OW; (cr, Ej,,) 


Js OE}, 


2 
= 2GE;,, + (x- c) (trE},.) I. (5.63) 
Furthermore, Equations (5.42), (5.55) and (5.60) yields the corresponding 
Cauchy stress 


1 
T = RMR’. (5.64) 


5.2.5 Field equations 


The chemical potential considering the coupling between diffusion and me- 
chanics is given by 
Ó n 


HrlorF)= = Me? HUR He (5.65) 


where L7"? and usd have been expressed in Equations (5.9) and (5.23), re- 
spectively. It should be noticed that, c and r in Equations (5.9) and (5.23) are 
replaced by cr and R, respectively, indicating these defined in the reference 


configuration. 
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For the coupling chemical potential i7", we consider two forms of it for the 


two different finite deformation elasticities: 


ETE RR 
Hg” (cr 1 m cm Div (a) 
=0 
ow? 
= OCR 


oye oss i gy. OF? 
oJ OCR oE! ` OCR 


1 oE! 
= Q(-.E:C:E “C: F°: 
FN C: E^) -/'C 3 


CR 
li: gu 
mo :E°) — Es iI 
1 1 
= Q(;E':C:E)- ;QuM', (5.66) 
l ayer ovy 
cp Fe = R _Dj R 
H, (cr, jog) OCR p (sei) 
= 
=0 
NE 
^. deg 
| ay? ays Ow; OE, 
i 0, Jes OE jog dcr 
1 QE; 
= D log :C: Ejog) + qe Flog ' ane 
1 Owe 
= 265 log ` zC: EU 37 EL, 
1 1 
= 265 Ejog: C: E Lu Qn 


9E: , 
Here, the term multiplied by 2 = or Ton represents the concentration depen- 


dence of the strain tensor in the chemical potential. As shown in the above 
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derivation, this term is equivalent to the external microforce term related to the 
Mandel stress in Anand's model [91, 94]. 
Combining Equations (4.9), (4.10), and (5.65) yields the coupled Cahn-Hilliard 


diffusion equation 
êr = Div (M(cn) Grad(ug"? + ut^ + eP) ) , (5.68) 


where M (c) has been defined in Equation (4.12). Again, replacing the term 
na by the term uuP"@lty, the coupled NSC diffusion equation with finite defor- 
mation elasticity can be obtained. 

Finally, based on the balance of momentum, as shown in Equation (3.21), the 


mechanical equilibrium condition is 


> 


div Tg — 6, (5.69) 


where the body force is neglected. 


5.2.6 Mathematical formulations of the spherically 
symmetric boundary value problem 


The coupled diffusion equation with elastic Green strain 


For the spherically symmetric boundary value problem, as mentioned in section 
5.1.4, only the radial component of displacement field is non-zero. As a result, 
the displacement gradient H expressed by Equation (3.3) can be rewritten in a 


diagonal form, only containing two different components Hj; and H22 = H33 


RE 
S 
5 

oO 


H = Gradüf(R,t) = (5.70) 


e e 
= 
= 


(en. 69.69] 
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Again, entities with subscript R indicate those defined in the reference config- 


uration. According to Equation (3.4), the deformation gradient F is expressed 


as 
ARAI 0 0 
F=H+I= 0 2-41 0 1 (5.71) 
0 0 m 


IZ 3731 


Combining Equations (5.37), (5.38), and (5.71) yields the elastic deformation 


gradient 


| | 9*1 0 0 
Diar EET 0 e+ 0 1 (5.72) 
0 D. I 


(en. 60.69] 


where f? = 4/1-- QO(cg — co) (see Equation (5.38)). Accordingly, based on 
Equations (5.43) and (5.72), the elastic Green strain tensor is expressed as 


E = 5 (FCF -1) 
Our 2 
a” = 0.5 0 0 
= 0 ED o5 0 
_ E 0 
+n? 
0 0 RHD gy 
2p {ér.€5 29} 


(5.73) 
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Substituting Equation (5.73) into Equation (5.52), the coupling energy density 


in terms of elastic Green strain tensor is given by 


V; (ca E^) 


; v 
GJ' |E^ : E 
| = 1—2v 


we] 


1 
= G(1+ (cr —co)Q) 7 7 1 


4( (e tp" i) 
2 UHR): 


ug \2 2 
(1+ 5) ) (5.74) 


Using Equation (5.55), the relationship between the second Piola-Kirchhoff 


stress tensor T° and the Cauchy stress tensor T can be rewritten as 


1 


T= FTF". (5.75) 


Combining Equations (3.12), (5.57), and (5.75), the first Piola-Kirchhoff stress 


tensor is given by 


Tr - I pepepel pT 
Je 
= pri 
2 
= 2GB'R (E+ («- 56) ee). (5.76) 
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As shown in Equation (5.76), we see the major diagonal entries of Tr are 
nonzero, and Tr,, = Tr,, due to the spherical symmetry. The components of it 


are therefore given by 


ı (ou 
TR, = G( (en - e)! (S +1) 
2 
ou 
Gr) 
(1+ (cr — c0)Q) 
2 
Ou, 
tae ($8 +1) (28 +1)? , 
1—2v | tler) Ale! i 
(5.77) 
l/u 
NS G(1-- (ca —09)9)* (+1) 
2 
ðu 
(3&1 
j 1 
(1 + (cr — c0)Q)* 
2 
Ou, 
m. (3g +1) | Q3 +1)” . 
1—2v \ Alert (1+ (cr-c0)Q)$ 
(5.78) 


Furthermore, due to spherical symmetry, the mechanical equilibrium (5.69) is 


expressed as 


Tr, 2 
ap Tg (TR Ten) =0. (5.79) 
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Now we derive the coupled diffusion equation of the spherically symmetric 
boundary value problem. First, the coupled chemical potential uf? (cg, E^), as 


shown in Equation (5.66), is expressed as 


2 
Our 
1 (1+2) 
uP(cr E) = GQ 14 t“, 
(1+ (er — co)Q)5 
2 2 
1+7% 
4l 1+ +) ;) 
(1+ (er c0)Q)* 
Ou, 
v 1 (1+ 98) 
+ 1+ 
1—2v 2 (1+(cr—co)Q)3 


Cai ) 


20 (14° (^ | (ee): 
1+ 


3 (1 + (cr — co)Q) (cr — c9)Q) 
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2 
-o (13g) 20 (14- #8)? 
360. Ares 
2, 
(+) 


(1+ (eg — c9)Q) 


ld 


(5.80) 


Substituting Equations (5.9), (5.23), and (5.80) into Equation (5.68), we can 


derive the dimensionless coupled diffusion equation 


0 = 


-20C o |ð 
OCR 
I or aR E 


(Ete a” ug! ges E). 


(5.81) 


For the detailed derivation of the final dimensionless coupled diffusion equa- 
tion, please see Appendix A.3.1. 

The boundary and initial conditions for the spherically symmetric boundary 
value problem are the same as those of the previous Cahn-Hilliard diffusion 
with SDT, as sketched in Fig. 5.4. Again, entities with subscript R in Fig. 5.4 


indicate those defined in the reference configuration. 
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3 
> m 

N | VALL: en = 

ur=0 Ta ig - 0 

—— «—— dc MES 

jg 7! Jg fig = const. 

i N Her M Gradeg ñR =0 
j| i 


Figure 5.4: Boundary and symmetry conditions for the spherically symmetric boundary value 
problem in the reference configuration [89]. 


The coupled diffusion equation with logarithmic elastic strain 


According to Equations (5.41) and (5.72), the right stretch tensor is expressed 
as 


U = vr'r 


1 
A1. 0 0 
1 
= (I-c(en—c0)Q) 5 0 R+l 0 
0 0 RH 


{ER 20 čo} 


(5.82) 


Therefore, combining the polar decomposition F^ = R°U* and the above Equa- 


tion (5.82) yields the rotation tensor 


R° =I. (5.83) 
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Furthermore, the spectral decomposition of U* can be obtained as 


3 
U = Y XMr&erfe. 
a=1 
= AE, G8, - A (1— 8, DER), (5.84) 
where 
Af = (1 oi tya 5.85 
ı = (1+(cr—co)Q) op I (5.85) 
1 
Yo = (1+ (cx 6)9) 5 C 1. (5.86) 


Thus, the logarithmic elastic strain tensor is expressed as 


m — In U* 
3 
377057 

a=1 

Inge Caf +1) : l 

i ; Ind. 4 +1) : 

A (UR 

Ü 0 In( gs C$ sb ) (en. 6o. 6o] 
(5.87) 


Here, the spatial logarithmic elastic strain tensor Ey, gH is equal to the logarith- 


mic elastic strain tensor EV, due to R^ = I. Substituting Equation (5.87) into 


log 
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Equation (5.53), the coupling energy density in terms of logarithmic elastic 


strain tensor is given by 


Verne) = GP Eti Eintr GE) 


ur +] 2 
= G(l-+(er-c0)Q) («( -)) 
(1+ (cr — co) Q)3 
2(n( Eu )) 
(1 + (cg — co)Q)? 
ve” (1( Sur +] ) 
1-2v (1+ (cr — c0)Q)3 


uR 2 
v2in( a JJ (5.88) 
(1+ (cr — co)Q)° 


Combining Equations (3.12), (5.64), (5.63) and (5.83), the first Piola-Kirchhoff 


stress tensor is given by 


J 
Tr = —-ME' 
Je 


2 
7 (2681, + (x- 50) (trE;,,) 1) FT. (5.89) 


The components of T are therefore given by 


Im = (terzo (ME +1) 


u an 
2GIn l 
EOT, 
2 OUR +1 
—=G){In OR 
is 52) ( Ga 


e ( 
(1+ (cg — c9)Q)3 


, (5.90) 
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Te, = (1+ (ence) o) (Aer) 


sew ( run J 

(1+ (cr—co)Q)3 
2 uR +] ) 
K--G In I 
+( A )( DET 


u gr )) 
(1 + (eg — c9)Q)3 


Now we derive the coupled diffusion equation of the spherically symmetric 
boundary value problem. First, the coupled chemical potential pP (cr, Ej, ). 


as shown in Equation (5.67), is expressed as 


(5.91) 
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(5.92) 


Substituting Equations (5.9), (5.23) and (5.92) into Equation (5.68), we can 
derive the dimensionless coupled diffusion equation 


0 = 


— OC d o zd 
29°R 2 mwp ed cp e 
R JF JR E ( R M(cR) (m + s tH, 5] : 


(5.93) 
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For the detailed derivation of the final dimensionless coupled diffusion equa- 


tion, please see Appendix A.3.2. 
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6 Phase-field modeling of 
sodium-ion battery particles 


In chapter 4, we introduced the phase field theory related to species insertion 
into or extraction from a cathodic particle where phase segregation possibly 
occurs during discharging or charging. In order to take the mechanical defor- 
mation into account, different elastic strain energy functionals were motivated 
from different mechanical theories in chapter 5. In this chapter, a phase-field 
model for NaFPO is studied for the first time (see also Zhang and Kamlah [89]). 
We implemented the model in COMSOL Multiphysics® for a spherically sym- 
metric problem of sodium insertion into or extraction from a cathodic particle 
made NaFPO. We first describe the determination of the material parameters 
for NaFPO. The model captures the important feature of phase segregation into 
a sodium-poor phase FePO, and a sodium-rich phase Naz/3FePO4. For this 
material, there is a visible difference for the concentration and stress between 
SDT and the finite deformation theories. Furthermore, we compare the two 
cathode materials NaFPO and LiFPO of lithium-ion batteries to each other in 
terms of phase changes and stresses, and show that, for the same stage of aver- 
age concentration, although the miscibility gap of NaFPO is smaller than that 
of LiFPO, the stresses in the cathode material NaFPO are higher in the phase 
segregated state. As a result, the suppression of phase segregation by the elastic 


strain energy is more easily achieved in NaFPO compared to LiFPO. 
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6.1 Material parameters of Na,FePO, 


We consider a typical active nanoparticle of size Ro = 150 nm for the particle 
radius. According to Nakayama et al. [120], the diffusion coefficient of Na- 
ions in NaFPO is an order of magnitude lower than that of Li-ions in LiFPO, 
which also can be confirmed by comparison to Zhu et al.[36]. The diffusion 
coefficient of LiFPO is 1 x 10714 m?/s [119, 121], therefore, we use the value 
Do = 1 x 107? m? /s as the diffusion coefficient of NaFPO. In the following, 
we will specify the values of the other parameters for NaFPO including Cmax, 
0, Oo, A. Q, and Ep. 


6.1.1 Determination of Cmax 


The maximum concentration Cmax with units of mol / m? means the maximum 
content of a certain species the host material can accept, which can be ex- 
pressed by [116] 


1 


E 6.1 
VoNa’ (6.1) 


Cmax = 


where Vo is the volume occupied by one species atom. This can be obtained by 


Vo = Veen (6.2) 
n 

with n being the number of atoms of the species that the host material can 
accept in a unit cell of volume Veel- 
According to Ref. [122], the structures of olivine LiFPO and olivine NaFPO 
are the same. Based on the crystal structure of LiFPO as shown in Ref. [123] 
and the lattice parameters in Table 6.1, we can calculate the maximum lithium 
concentration of LiFPO by Equation (6.1), which gives 2.29 x 10^ mol/m?, 
matching the value used by Zeng and Bazant [119], and Delacourt and Safari 


[124]. Therefore, in the same way, the maximum sodium concentration of 
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6.1 Material parameters of Na,FePO4 


NaFPO can be obtained as Cmax = 2.1 x 10^ mol /m?, which is smaller than 


that of LiFPO due to the larger unit cell volume V..;; for NaFPO. 


LiFPO NaFPO 
Crystal structure Olivine Olivine 
a (A) 10.332 [125] 10.4063 [122] 
b(A) 6.010 [125] 6.2187 [122] 
c (A) 4.692 [125] 4.9469 [122] 
Veell " 291.4 [125] 320.13 [122] 
n 4 4 


Table 6.1: Lattice parameters of LiFPO and NaFPO. 


6.1.2 Determination of o, and a 


According to Lu et al. [37], at room temperature, the phase diagram of olivine 
NaFPO consists of two regions. For 0 « x « 2/3, phase segregation into a 
sodium-poor phase FePO4 and a sodium-rich phase Na» /3FePO4 is found to 
be favorable, which means that this is a two-phase region. (Following common 
praxis, we replace in chemical formulas dimensionless concentration c by x.) 
Based on the unit cell parameters and volume as functions of x in NaFPO as 
shown in Ref. [37], two different values for the unit cell parameters and vol- 
ume, respectively, can be identified when x increases to around 0.08, which is 
the manifestation of the occurrence of the phase segregation. Therefore, we can 
obtain the normalized sodium concentration for the initiation of phase segre- 
gation, which is around 0.08. In contrast to LiFPO, where transformation from 
FePO4 into LiFePO, occurs directly, the system of NaFPO goes through an 
intermediate state at Na? 3 FePOA4. For 2/3 < x < | there is the solid-solution 
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phase Na,FePO, which is a single-phase region. With regard to the homoge- 


neous free energy density in the two-phase region of NaFPO, we assume 


Q 
w""P(c) — kpTserNACmax CE 27 


I) 


This is the classical homogeneous free energy density function for a two phase 


(6.3) 


material [27, 85], which we have formulated such that it is limited to the range 
0 « x « 2/3. Consequently, we fit the homogeneous free energy density (6.3) 
to the above experimental data [37] with respect to the parameters a and 06, 
as shown in Fig. 6.1. The best fit was achieved with o = 5 and o5 = —15 at 


room temperature. 
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Figure 6.1: Fit of the normalized homogeneous free energy density in the two-phase region to the 
experimental data from Ref. [37] [89]. 
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6.1 Material parameters of Na,FePO4 


Fig. 6.1 shows the fitting result with the normalized homogeneous free energy 
density in the two-phase region yw”? plotted versus normalized concentration 
C. It can be seen that i"? exhibits a doublewell structure, such that two dif- 
ferent relative minima A and B occur, characterizing two phases of different 
solubility limits cq and cg respectively, of the sodium concentration. The two 
minima A and B correspond to the two phases FePO4 and Na» j3FePOA, re- 
spectively. The two zones AC and BD are the "nucleation zones", and phase 
segregation occurs upon sufficient disturbance of the system's equilibrium. In 
the inner zone of concavity between the two points of inflection C and D, 
which is denoted the “spinodal decomposition zone", homogeneous sodium 
concentration states are unstable and phase segregation is initiated in any case. 
The Maxwell construction represents the volume fractions of a phase segre- 
gated system. The inflection points correspond to the second order derivative 
of y""? with respect to c, as shown in Fig. 6.2, and the inflection point C 


matches the experimental data from Ref. [37]. 
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Figure 6.2: Second order derivative of y"? with respect to c [89]. 
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6.1.3 Determination of A 


The so-called reference value Ag for A can be estimated as [116] 


2 
3 


= ica 
ym e uc (wiz) 


3 ; (6.4) 


and the relation between Ao and A is 
À = «Ào, (6.5) 


where & is constant. 

Using Equation (6.4), we obtain Ag = 9.1 x 107? m? for NaFPO. In absence 
of any further information, we pick œ = 20 as the value for NaFPO which is 
the one digit number in the order of magnitude of the values for LiFPO and 
LMO, see Table 6.2. Therefore, we get A = 1.8 x 10-17 m? for NaFPO. The 


values of À for three cathode materials are shown in Table 6.2, as well. 


Ào a A 
LMO |3x10-? (m?) 23 7x 10-18 (m?) [ ua 
LiFPO | 5.2 x 10-19 (m?) 17 8.8 x 10-18 (m?) [5 
NaFPO | 9.1 x 107? (m?) 20 1.8 x 107 (m?) 


Table 6.2: Gradient energy coefficients of three cathode materials. 
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6.1.4 Determination of Q 


The partial molar volume © plays a role analogous to a thermal expansion 
coefficient, meaning we can calculate Q by the relation &° = 1/3Q(c — co)I. 
Thus, based on Table 6.3, the partial molar volume of LiFPO is obtained as 
0.022 x 3 6 3 
Q = — — =2.9x 10” m /mol, (6.6) 
(1 = 0) Cmax / 
which matches the value used by Song et al. [117]. Therefore, in the same 


way, the partial molar volume of NaFPO can be calculated as 


0.041 x 3 


————— - 8.8 x 1075 m? /mol. (6.7) 
(3 - 0) Cmax 


This value of © is larger than that of LiFPO, which is consistent with the fact 
that sodium has a larger cation radius than lithium as shown in Table 1.1. It 


should be noticed that we do not consider the possibility of a concentration 


dependence of Q. 
Li,FePO4 Nay,FePO4 
x 0cx«l 0<x< 2/3 
volume expansion 6.8% [37] 12.8% [37] 
strain 0.022 0.041 


Table 6.3: Volume change of LiFPO and NaFPO. 


6.1.5 Determination of Eo 


Young's modulus of the host material FePO4 is 125 GPa, Young's modulus of 
LiFePO, is 123.9 GPa, and Young's modulus of the cathode material LiFPO is 
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124.5 GPa, which is the average of the values for FePO4 and LiFePO4 [126]. 
Based on the above data, due to the larger ionic radius of the Na-ion compared 
to that of a Li-ion, we regard the two digit estimate 120 GPa as Young's mod- 
ulus of NaFPO. The material parameters for the two cathode materials LiFPO 
and NaFPO are summarized in Table 6.4. 


Parameter LiFPO NaFPO 

0 4.5 [50] 5 

0 -9 [50] -15 

À 8.8 x 1071? (m?) 1.8 x 107 (m?) 

Do 1 x 10714 (m? /s) 1 x 10-15 (m?/s) 
Cmix 2.29 x 10^ (mol /m?) 2.1 x 104 (mol /m?) 
Q 2.9 x 1076 (n? /mol) 8.8 x 1075 (m? /mol) 
Eo 124.5 (GPa) 120 (GPa) 

v 0.25 [126] 0.25 [126] 


Table 6.4: The material parameters for the cathode materials LiFPO and NaFPO. 


6.2 Numerical implementation 


The resulting set of equations has been implemented in COMSOL 
Multiphysics? for solution by the finite element method. For the numerical 
methods, since the Cahn-Hilliard equation is a fourth order partial differential 
equation for the species concentration field, the standard finite element method 
with C°-continuous Lagrange basis functions is not sufficient for discretiza- 
tion. There are several methods to overcome these numerical difficulties. For a 
straightforward treatment of the high-order operator, Zhao et al. [97] employed 


the isogeometric analysis, which allows for the employment of a wide range of 
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smooth, higher-order basis functions, providing global C!-continuity. An al- 
ternative approach is the use of split-methods or mixed-methods to reduce the 
fourth order equation into two second-order equations. For example, Miehe et 
al. [127] considered the chemical potential as an additional degree of freedom 
to derive the mixed form for C°-continuous basis functions. Huttin and Kamlah 
[27] introduced the second derivative of the species concentration as unknown 
function to split the Cahn-Hilliard equation into two second-order equations, 
avoiding the problem with the computation of the logarithmic entropic terms 
in the homogeneous free energy density as shown in Equation (4.4). In ad- 
dition, nonlocal species concentration [128] or micromorphic species concen- 
tration [94] is introduced as an extra degree of freedom to obtain an efficient 
implementation with C°-continuous Lagrange basis functions. In our imple- 
mentations, we use the second derivatives of the species concentration as an 


additional degree of freedom. 


6.3 Cahn-Hilliard model without 
mechanics 


We consider the quasistatic insertion and extraction of sodium for a particle of 
NaFPO at C — 0.001, see the boundary condition (5.30). With this C-rate, we 
study the behavior for dynamic, i.e. continuous insertion very close to a se- 
quence of equilibrium states. In this way, the system is allowed to move along 
a path of relaxed quasiequilibrium states. In our simulations, we exclusively 
focus on the two-phase region of NaFPO (0 « x « 2/3). The average concen- 
tration cay, also called “state of charge” (SOC) is defined as cayg = fg cdV /V. 
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Figure 6.3: Normalized average system free energy Y ug and normalized homogeneous free energy 
density U""? as function of Cavg and c, respectively [89]. 


In Fig. 6.3, the markers represent values of the normalized average system free 


energy 


IL I nw Trga = 
Yan = 3 n "P + 81 (grad c))dV (6.8) 


during insertion and extraction, respectively, as function of Cavg. For demon- 
stration purposes, the plot of the normalized homogeneous free energy density 
vs. dimensionless concentration is entered, as well. Markers on the normal- 
ized homogeneous free energy density curve correspond to homogeneous states 
whereas markers nearby the path of the Maxwell construction correspond to 


phase segregated states, as they are illustrated in Figs. 6.4 and 6.5. Fig. 6.4 
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Figure 6.4: Normalized sodium concentration © versus normalized radial coordinate r/Ro for dif- 
ferent time instants during sodium insertion [89]. 


shows the normalized sodium concentration along the radial coordinate at dif- 
ferent time instants during sodium insertion. We find that at the beginning of 
sodium insertion the concentration is homogeneous and the corresponding con- 
centration is 0. Once Cavg gets close to 0.08, which corresponds to the inflection 
point C in Fig. 6.1, phase segregation is initiated. Two different phases can be 
recognized in the particle, namely the sodium-poor phase FePO, in the center 
corresponding to the first minimum A and the sodium-rich phase Na» /5FePO4 
at the outside corresponding to the second minimum B. A smooth but very 
narrow interface with concentration changing rapidly but continuously sepa- 
rates them. When Cavg grows up to 2/3, the intermediate phase Na»;5FePO4 


occupies all of the particle. 
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ferent time instants during sodium extraction [89]. 


Next, we consider the sodium extraction case shown in Fig. 6.5. At the begin- 
ning the system is in a homogeneous state and the corresponding concentration 
is 2/3. Once Cayg is reduced to around 0.6, which corresponds to the inflection 
point D in Fig. 6.1, phase segregation is initiated during sodium extraction. 
In contrast to the insertion case, the sodium-poor phase FePO, is at the out- 
side while the sodium-rich phase Na» 5FePO is in the center. At the end of 
sodium extraction, the sodium-rich phase Na» /FePO, vanishes. It should be 
noticed that although to the particle surface a mass flux is applied, the concen- 
tration on the particle surface nearly stays at a constant value. This is due to 
the fact that because of the extremely low C-rate of C — 0.001 for quasistatic 
insertion and extraction of sodium, the system moves along a path of relaxed 


quasiequilibrium states. 
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6.4 Cahn-Hilliard model with mechanics 


Concerning mechanics, we first study the coupling between the Cahn-Hilliard 
theory and SDT, and then compare the different small and finite deformation 
theories. In order to study the coupling effect through the coupling energy on 
sodium diffusion and stress in the particle, we introduce a normalized Young's 


modulus as 


E- a (6.9) 

TE : 
where Eo is the value of Young's modulus shown in Table 6.4. In the figures, c, 
cr denote concentration of small and finite deformation elasticity, respectively, 


and Ty = 1/3T,;; is the hydrostatic stress in terms of the Cauchy stress. 


6.4.1 Cahn-Hilliard model with small deformation 
theory 


Figs. 6.6 and 6.7 show normalized sodium concentration and hydrostatic stress 
along the radial coordinate at a stage of average concentration of Cavg = 0.5 
during sodium insertion based on SDT for different E. We find that the differ- 
ence in concentration between the two phases is reduced as E increases, which 
means that the solid solution limits of FePO4 and Na, FePO, are gradually ex- 
tended into the range of phase segregated states. In other words the coupling 
effect reduces the miscibility gap. Correspondingly, the hydrostatic stress mag- 
nitudes increase as E increases. The hydrostatic stress is constant and tensile in 
the low concentration phase, and it drops rapidly across the interface to become 
compressive and constant in the high concentration phase. This change of sign 
is a result of mechanical equilibrium between the inner core and the outer shell 
of the particle. Note that on the other side for E = 1, 
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the system is homogeneous and stress-free. Even though the average concen- 


tration value lies in the spinodal decomposition zone of W”""P, phase segrega- 


tion does not occur. 


NaFPO, R,=150 nm 
C -0.001, c,,,=0.5 
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Figure 6.6: Normalized sodium concentration © versus normalized radial coordinate r/Ro at Cavg = 
0.5 based on SDT for different E during sodium insertion [89]. 
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Figure 6.7: Hydrostatic stress Tg versus normalized radial coordinate r/Ro at Cavg = 0.5 based on 
SDT for different E during sodium insertion [89]. 


As illustrated in Fig. 6.8, for a high value of normalized Young's modulus, 


the cost of the system coupling energy V? at a phase segregated state would 
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be too high so that as a result of the competition between the different contri- 
butions to the system free energy the system stays in the homogeneous state 
in order to minimize the total system free energy. Therefore, the presence of 
a coupling energy in the total system free energy can lead to suppression of 
phase segregation. For the detailed discussion of the influence of mechanics 


on the miscibility gap, see chapter 7. 
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Figure 6.8: Normalized system coupling energy W^? versus average concentration Cave based on 
SDT for different E during sodium insertion. Here Y? = fy ?dV [89]. 


6.4.2 Comparison of different mechanics theories 


Figs. 6.9 and 6.10 show the comparison of normalized sodium concentration 
and hydrostatic stress at a stage of average concentration of Cayg = 0.5 dur- 


ing sodium insertion for different mechanics theories with a value of E = 0.3: 
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SDT, elasticity based on elastic Green strain, and elasticity based on logarith- 
mic elastic strain. According to the experimental work of Xiang et al. [39], 
driven by the coupling energy, the sodium-rich phase decreases its sodium con- 
tent to about Nao.sFePO4, while the sodium-poor phase increases its sodium 
content to about Nao.osFePO4. These experimental results [39] are shown as 
dashed horizontal lines in Fig. 6.9. It can be recognized that there is a visible 
difference for the sodium concentration plots of SDT and the finite deformation 
elasticity formulations, with the difference between the two phases being less 
pronounced for SDT, while the sodium concentration plots of two finite defor- 
mation theories are almost the same. Compared to the experimental result [39], 
we can find that the sodium concentration in the high concentration phase ob- 
tained from two finite deformation theories matches the experimental value for 
Nao.«FePOs. but there is a clear deviation from the result by SDT. On the other 
hand, the sodium concentration in the low concentration phase obtained from 
all three mechanics theories is always slightly below the experimental value 
for Nao.ogFePO.4. In Fig. 6.10, the hydrostatic stress plots are clearly different 
for SDT on the one side and the two finite deformation theories on the other. 
There are higher stresses in the two phases for SDT. This is attributed to the 
fact that SDT does not account for the volume swelling of the particle during 
sodium insertion. Comparing the two finite deformation theories to each other, 


it is found that there is a negligible difference in stress between them. 
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Figure 6.9: Normalized sodium concentration versus normalized radial coordinate at Cayg = 0.5 
with E — 0.3 during sodium insertion for different mechanics theories. The normal- 
ized radial coordinates r/Ro and R/Ry represent SDT and finite deformation theory, 
respectively. The dashed horizontal lines represent the experimental results from Ref. 


[39] [89]. 


As mentioned before, the contribution of the coupling energy to the total sys- 
tem free energy can lead to suppression of phase segregation. Here, we discuss 
the critical value of normalized Young's modulus E, of NaFPO above which 
phase segregation can not arise in the particle. Fig. 6.11 shows the normalized 
sodium concentration along the radial coordinate for various values of E dur- 
ing sodium insertion based on different mechanics theories. For SDT, it can 
be seen that phase segregation can arise at a stage of average concentration of 
Cavg = 0.33 when E = 0.385 but sodium concentration is homogeneous over 
the particle when E = 0.386, so E, is 0.385 for SDT. In the same way, we find 
E. = 0.409 for the two finite deformation elasticity formulations, which is a 
little larger than that for SDT. 
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Figure 6.10: Hydrostatic stress versus normalized radial coordinate at Cavg = 0.5 with E=0.3 
during sodium insertion for different mechanics theories [89]. 


6.5 Comparison with the cathode material 
Li,FePO, 


We now compare the two cathode materials NaFPO and LiFPO in terms of 
phase changes and hydrostatic stress. Figs. 6.12 and 6.13 show normalized 
concentration and hydrostatic stress along the radial coordinate for various 
particle radii Ro for the two cathode materials during insertion based on dif- 
ferent finite deformation theories. Here, we consider two stages of average 
concentration, one is Cayg = (Ca -- cg)/2 which is the center of the combined 


nucleation and spinodal zone, the other is Cavg = 0.5. It should be noticed that 
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Figure 6.11: Normalized sodium concentration versus normalized radial coordinate for various E 
during sodium insertion based on different mechanics theories [89]. 


(Ca +cp)/2 of LiFPO is just equal to 0.5, and which of NaFPO is 0.333. In 
Fig. 6.12, we recognize that there is an obviously reduced miscibility gap for 
a NaFPO particle compared to a LiFPO particle. This is consistent with the 
experimental fact that NaFPO has a two-phase region in the range 0 « x « 2/3 
while phase segregation for LiFPO occurs in the range 0 < x « 1 [37]. For 
Cave = (Ca teg) /2, the two cathode materials show almost the same position 
of the center of the interface region, which is located close to the particle sur- 
face due to the spherical symmetry. Additionally, for Cavg = 0.5, the interface 
location in a NaFPO particle is more close to the particle center, and, further- 


more, the concentration in the two phases is the same for different particle radii. 
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As shown in Fig. 6.13, for ca; = (Ca + cg)/2, both the tensile stresses in 
the inner core and the compressive stress magnitudes in the outer shell of a 
NaFPO particle are larger than those in a LiFPO particle. When Cayg increases 
from 0.333 to 0.5, the tensile stresses in the inner core of a NaFPO particle 
become larger, while the compressive stress magnitudes in the outer shell of a 
NaFPO particle decrease, which is even smaller that those in a LiFPO particle. 
Indeed, the particle shell made of high-concentration phase becomes thicker 
when Cavg shifts from 0.333 to 0.5 (please see Fig. 6.12), so the particle core 
made of low-concentration phase is in mechanical equilibrium with a thicker 
shell made of high-concentration phase. As a result, both the tensile stresses 
in the inner core and the compressive stresses in the outer shell of a NaFPO 
particle increase. Therefore, for the same stage of average concentration, al- 
though the miscibility gap of NaFPO is much smaller than that of LiFPO, the 
stresses in a NaFPO particle are larger compared to those in a LiFPO particle. 
This is due to a larger expansion during phase changes for NaFPO as shown in 
Fig. 6.14 describing the plots of volume ratio J at the particle surface. Here, 
it should be mentioned again that our simulations only consider the two-phase 
region of NaFPO (0 < x < 2/3). We can find that the volume expansion from 
FePO4 to Naz/3FePOg is quite large, namely about 12.8%, which is nearly 2 
times that for LiFPO changing from FePO4 to LiFePO4. This is consistent 
with the reports [37-39]. Correspondingly, as shown in Fig. 6.15, the maxi- 
mum hydrostatic stress magnitude in a NaFPO particle during phase changes 
is always larger than that in a LiFPO particle. The larger stresses in a NaFPO 
particle may explain the existence of a wide range of solid solution Na,FePO4 
(2/3 « x « 1) to avoid even higher stresses in the NaFPO particle. Casas- 
Cabanas et al. [38] conclude that the larger stresses in a NaFPO particle can be 
the explanation for the existence of an intermediate phase. Indeed, the forma- 
tion of an intermediate phase acts as a buffer between FePO4 and NaFePO4 
providing elasticity to the structure. On the other hand, the critical value of 
normalized Young's modulus E; of LiFPO for surpression of phase segrega- 


tion is equal to 1 (please see Fig. 7.5 in chapter 7), which is larger than that of 
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NaFPO, as shown before. Therefore, it 1s easier for NaFPO to reach a massive 
coupling energy to suppress phase segregation compared to LiFPO due to the 
larger volume change of a NaFPO particle. 

As can be seen in Fig. 6.15, the LiFPO particle surface expands more seriously 
around Cavg = 0.15. This is due to phase segregation being already initiated 
for LiFPO but not yet for NaFPO which is still in the homogeneous state. As 
a result, the maximum hydrostatic stress magnitude in a LiFPO particle is far 
greater than that in a NaFPO particle at this insertion state. At c;,4 = 2/3, 
the maximum hydrostatic stress magnitude in a NaFPO particle is close to zero 
while the LiFPO particle still displays high stresses. This is attributed to the in- 
termediate phase Na» /FePO, occupying then the whole particle while LiFPO 


is still in a phase segregated state. 
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Figure 6.12: Normalized concentration Cr versus normalized radial coordinate R/Ro with E = 0.3 
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for various Ro during insertion for two cathode materials based on different finite 
deformation elasticity formulations. Two stages of average concentration are consid- 
ered: Cayg = (Ca +cg)/2 and cay; = 0.5, where (ca +cg)/2 of LiFPO is just equal to 
0.5, and which of NaFPO is 0.333. 


6.5 Comparison with the cathode material Li,FePO4 
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Figure 6.13: Hydrostatic stress Ty versus normalized radial coordinate R/Ro with E = 0.3 for 
various Ro during insertion for two cathode materials based on different finite defor- 
mation elasticity formulations. Two stages of average concentration are considered: 
Cave = (Ca +cp)/2 and Cayg = 0.5, where (Ca + cg)/2 of LiFPO is just equal to 0.5, 
and which of NaFPO is 0.333. 
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Figure 6.14: Volume ratio J at the particle surface versus ca,g with E = 0.3 for various Ro during 
insertion for two cathode materials based on different finite deformation elasticity 
formulations [89]. 
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Figure 6.15: Maximum hydrostatic stress magnitude |Tg max| versus Cayg with E — 0.3 for various 
Ro during insertion for two cathode materials based on different finite deformation 
elasticity formulations [89]. 
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7 Particle size and average 
concentration dependent 
miscibility gap 


In the previous chapter 6, we just focus on the bulk particle size of the cath- 
ode material. As mentioned before, the thermodynamics of phase segregation 
in nanoscale particles is distinctly different from bulk materials, and the parti- 
cle size miscibility gap plays a nontrivial role in the performance of nanoscale 
insertion materials. In contrast to our common thermodynamic knowledge, 
the miscibility gap in nanoscale insertion material also depends on the average 
concentration. In this chapter, the particle size and average concentration de- 
pendent miscibility gap of the nanoscale insertion materials LMO, LiFPO, and 
NaPPO are investigated during insertion, using a coupled phase-field model 
based on the Cahn-Hilliard theory and finite deformation elasticity (see also 
Zhang and Kamlah [96]). For the mechanical part, two finite deformation elas- 
ticity formulations, based on elastic Green strain and logarithmic elastic strain, 
respectively, have been compared in the previous chapter 6. We found that, al- 
though the results from the two formulations of Hooke's law with the geometric 
nonlinearity of finite strain are almost the same, according to our experience, 
elasticity based on logarithmic elastic strain is numerically more efficient than 
elasticity based on elastic Green strain. The stress strain relation for LiFPO 
is also described by elasticity based on logarithmic elastic strain in the work 
by Di Leo et al. [94]. Therefore, a logarithmic elastic strain law is employed 
in this chapter. We implemented the model in COMSOL Multiphysics® for 
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a spherically symmetric problem of species insertion into cathodic particles 
which are respectively made of LMO, LiFPO, and NaPPO. We first compare 
three bulk cathode materials in terms of phase changes and stresses, and then 


investigate their particle size and average concentration dependent miscibility 


gap. 


7.1 Material parameters 


In order to study the particle size and average concentration dependent misci- 
bility gap in different nanoinsertion cathode materials, the three cathode mate- 
rials LMO, LiFPO, and NaFPO are chosen in the simulations. The material pa- 
rameters for these three cathode materials are summarized in Table 7.1, where 
the material parameters of NaFPO have been determined in chapter 6. We vary 
the nanoparticle radius to study the influence of particle size on the miscibility 
gap. Depending on the temperature, LMO possibly exhibits phase segregation 
in the range of values 0 < x « 1 where the crystalline host structure remains 
cubic spinel [132-134]. Some theoretical calculations of LMO have been car- 
ried out based on an assumption that Li-ordering phase exists at x — 0.5 around 
room temperature [23, 135, 136], but this does not agree with the experimental 
literature that suggest no Li-ordering phase [137]. At x — 1, LMO undergoes 
a phase transition where the crystalline material becomes tetragonal due to the 
Jahn-Teller distorsion regarding the manganese ions [138]. For LMO, we fo- 
cus on the phase segregation that happens at room temperature in the range of 
values 0 « x « 1, neglecting the effect of lithium ordering. The cubic to tetrag- 
onal transition at x — 1 is also not taken into account since the state of charge 
is rarely driven up to x — 1 unless delivering a large current density [139]. 

At room temperature, the miscibility gap of bulk LiFPO encompasses nearly 
the entire lithium composition range of the material [50] compared to a rel- 
atively small miscibility gap of bulk LMO bounded by solid solutions states 
between Lig 12Mn5O4 and Lio sgMn5O4 [116]. However, bulk NaFPO shows 
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Parameter | LMO LiFPO NaFPO 

0 2.5 [116] 4.5 [50] 5 

Qh -5.2 [116] -9 [50] -15 

À 7 x 107'8(m?) | 8.8x 10775 (m?) | 1.8 x 1077 (m?) 
[116] [50] 

Do 7.08 x | 1x 1074 (m?/s) | 1x 107} (m?/s) 
107? (m? /s) [121] 
[129] 

Cage 2.29 x | 2.29 x 104 | 2.1 x 10^ 
10^ (mol /m?) (mol/m?) [124] | (mol/m?) 
[130] 

Q 3497 x 10.29 x 10$5|88 x 10% 
(m? /mol) [130] | (m3/mol) [117] | (m? /mol) 

Eo 93 (GPa) [131] 124.5 (GPa) | 120 (GPa) 

[126] 
v 0.3 [131] 0.25 [126] 0.25 [126] 


Table 7.1: The material parameters for three cathode materials. 


the smallest miscibility gap bounded by solid solutions states between FePO4 
and Na» j5FePO4 [37] among these three cathode materials. 

The three cathode materials show different volume changes during species in- 
sertion. For LMO, the volume expansion of Mn5O4 upon lithiation to LiMn204 
reaches about 7.396 [140]. According to Lu et al. [37], the volume expansion 
of FePO, upon sodiation to NaFePO; is quite large, namely about 17%, and 
even the volume expansion from FePO4 to Na»/5FePO, is still quite large 
(about 12.8%), which is nearly 2 times that for LiFPO (about 6.8%) changing 
from FePOy to LiFePO.. The previous chapter 6 reveals that the difference for 
the concentration and stress plots between SDT and finite deformation elas- 


ticity can not be neglected for NaFPO. Therefore, finite deformation elasticity 
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is selected in order to represent all the three cathode materials by the same 


deformation theory. 


7.2 Comparison of the three bulk cathode 
materials 


We will consider quasistatic insertion of species into nanoparticles for three 
cathode materials at C = 0.001. Again, Baye = f. By VRdVn / Vn is the normal- 
ized average system free energy, see also Equation (5.48). We compare three 
bulk cathode materials in terms of phase changes and stresses, and a typical 
active particle radius of Ro — 200 nm is considered. 

Fig. 7.1 shows the values of the normalized average system free energy for 
the three cathode materials during insertion by the markers, including the pure 
diffusion case without mechanical deformation and the mechanically coupled 
diffusion case, as function of c;,,. For demonstration purposes, the plots of the 
normalized homogeneous free energy density versus the normalized concentra- 
tion are also entered as line plots. The species-poor phase and the species-rich 
phase are characterized by the solubility limits cq and cg, respectively, of con- 
centration. Here, the solubility limits are the lower and upper bounds of the 
concentration range of phase segregated states. For example, in Fig. 7.1, the 
solubility limits can be identified where the normalized average system free 
energy states do not coincide with the plot of the homogenous multiwell po- 
tential. The difference in the value of W”"P at a stage of average concentration 
of Cave = (Ca +cg)/2, and the value of y"? at either of the solubility lim- 
its cy or cg represents the maximum energy barrier AU?"". As mentioned 
in section 6.1.2, the two ranges between the respective point of tangency de- 
fined by the Maxwell construction and the corresponding inflection point are 
the "nucleation zones", and the inner zone of concavity between the two points 


of inflection is called the “spinodal decomposition zone". Again, markers on 
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the normalized homogeneous free energy density curve correspond to homo- 
geneous states whereas markers nearby the path of the Maxwell construction 
correspond to phase segregated states. In particular, at a stage of average con- 
centration of Cgyg = 0.5, all three cathode materials show phase segregated 


state, as illustrated in Fig. 7.2 (a). 


C = 0.001, insertion 
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Figure 7.1: Normalized average system free energy VP iig and normalized homogeneous free energy 
density U/""? as function of Cayg and c, respectively, for the three cathode materials and 
aradius of Ro = 200 nm. The dashed straight lines represent Maxwell construction, and 
E — 0.3 is considered in the mechanically coupled diffusion case [96]. 


We first focus on the pure diffusion case. In Fig. 7.1, we find that phase 
segregation is initiated for NaFPO once c;,4 approaches the inflection point 
at Cayg = 0.08, which is earlier compared to LiFPO (c4,, = 0.13) and LMO 
(Cave = 0.26). Also, it can be seen that LiFPO has the largest “spinodal de- 


composition zone" among the three cathode materials. In Fig. 7.2 (a), we 
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also find that the miscibility gap, when neglecting mechanics, is the smallest 
for NaFPO, which is about 2/3 of that of LiFPO. However, if the contribution 
from the coupling energy is considered (E — 0.3), the miscibility gap of LMO 
is significantly reduced. We can see that in this case the related normalized 
average system free energy shown in Fig. 7.1 coincides with the respective 
normalized homogeneous free energy density even in most of the "spinodal 
decomposition zone". Hardly visible, phase segregation, i.e. a minor deviation 
from yw”? takes place in the neighborhood of cgyg ~ 0.5, only. The result- 
ing miscibility gap is even smaller than that of NaFPO, even if mechanics is 
accounted for the latter material. The physical explanation for this behavior 
of LMO is that the cost of the elastic strain energy at a phase segregated state 
would be too high so that as a result of the competition between the different 
contributions to the system free energy the system stays in the homogeneous 
state in order to minimize the total system free energy. The effect of mechan- 
ics, i.e. the system coupling free energy, is the weakest for LiFPO. Indeed, as 
shown in Fig. 7.1 (a), in contrast to the other two materials, when considering 
mechanics, the markers of LiFPO are still close to the path of the Maxwell 
construction, and the miscibility gap of LiFPO is just slightly reduced in Fig. 
7.2 (a). As a result of the larger miscibility gap, which by the partial molar 
volume translates into a larger strain mismatch, the tensile stresses in the inner 
core of a LiFPO particle are larger compared to those in a LMO particle, and 
the compressive stress magnitudes in the outer shell of a LiFPO particle are 
also larger than those in a LMO particle, as shown in Fig. 7.2 (b). On the other 
hand, as mentioned in section 6.5, although the miscibility gap of NaFPO is 
smaller than that of LiFPO, both the tensile stresses in the inner core and the 
compressive stresses in the outer shell of a NaFPO particle are larger than those 
in a LiFPO particle at c;,, = 0.5. This again is due to a larger expansion, i.e. 


partial molar volume, during phase segregation for NaFPO. 
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Figure 7.2: (a) Normalized concentration versus normalized radial coordinate at cgyg = 0.5 during 
insertion for three cathode materials of Ro = 200 nm. The solid lines represent the 
pure diffusion case and the dashed lines represent the mechanically coupled diffusion 
case of E = 0.3. (b) Hydrostatic stress Ty versus normalized radial coordinate R/Ro 
at Cayg = 0.5 during insertion for three cathode materials (Ro = 200 nm and E =0.3) 
[96]. 


7.3 Particle size dependent miscibility gap 


In order to investigate the impact of the particle size on the miscibility gap for 
the three cathode materials, we first exclude the effect of the mechanics. The 
solubility limits of the material are determined by taking the minimum and 
maximum concentrations for which phase segregated states occur [50]. 

Fig. 7.3 shows the solubility limits as a function of particle radius at a stage of 
average concentration of Cavg = (Ca + cg )/2 for three cathode materials. The 
motivation for focusing on this specific value of average concentration is based 
on observations in our simulation results, that the miscibility gap shrinks sym- 
metrically about the center of the combined nucleation and spinodal zone. 


Consequently, we assume that phase segregation is suppressed completely, 
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when it is suppressed for Cayg = (Ca +cp)/2. Here, cq and cg are the size 
independent solubility limits found for sufficiently large particles. 

We find that a reduction of the particle radius leads to a shrinking of the mis- 
cibility gap. In particular, for each material a specific critical particle radius 
is found below which phase segregation is inhibited. The critical radius of 
LMO (R. = 10.9 nm) is larger than those of 6 nm for LiFPO and 6.4 nm for 
NaFPO, respectively. The above predicted critical particle sizes match the fact 
that phase segregation would be suppressed when the particle size becomes 
comparable to the interface thickness. According to Cahn and Hilliard [29], 
the interface thickness d is proportional to the square root of the gradient coef- 
ficient A and the inverse of the maximum energy barrier , i.e d (A /AW""P) 2. 
Therefore, LMO has the largest critical particle size among the three cathode 
materials due to a very tiny AU"? shown in Fig. 7.1. It should be mentioned 
that our value of R. = 6 nm for LiFPO is a little bit smaller than the estimate 
value R. = 7.5 nm from Meethong et al. [45], but it is close to the critical 
radius (around 5 nm) obtained from a 3D comprehensive phase-field model ac- 
counting for facet dependent surface wetting [51]. On the other hand, different 
materials show different particle size dependent miscibility gap behavior. For 
LMO a reduced miscibility gap is found below a particle radius of 55.5 nm. 
The particle size dependent miscibility gap becomes more and more obvious 
as the particle size decreases, and strongly varying solubility limits are found 
below a particle radius of 35 nm. Compared to LMO, the threshold values 
of the particle radius for the size independent miscibility gap are smaller for 
LiFPO and NaFPO, which are 20.5 nm and 17.5 nm, respectively. Therefore, 


LMO has the broadest particle size range for varying solubility limits. 
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Figure 7.3: Solubility limits as a function of particle radius at cayg = (co + cg) /2 during insertion 
for three cathode materials [96]. 


With the help of Fig. 7.4, we take LiFPO as an example to explain the shrink- 
ing of the miscibility gap upon particle size reduction. As the particle size 
reduces, the interface region, the thickness of which is determined by the ma- 
terial constant À, covers a larger fraction of the material. Furthermore, due to 
the spherical symmetry, the position of the center of the interface region for 
Cave = (Ca +cp) /2 ~ 0.5 is located close to the particle surface for all values 
of radius. Consequently, below a radius of 20.5 nm the concentration value at 
the surface progressively is reduced from an initial value of cg. As a result 
the miscibility gap shrinks at the same average concentration to accommodate 
in a particle of reduced size the concentration gradient, the slope of which is 
basically given as a material constant. Due to the normalization of the radial 
coordinate in the form r/Ro in Fig. 7.4, this effect translates into an apparent 


decrease of the concentration gradient, resulting in the same decrease of the 


111 


7 Particle size and average concentration dependent miscibility gap 


surface concentration. Once the particles reduces to be below the critical parti- 
cle size, the cost of the gradient energy ‘P84 at a phase segregated state would 
be too high so that the system would stay in a homogeneous state in order to 
minimize the total system free energy due to the competition between the dif- 
ferent contributions to the system free energy, namely multiwell potential and 


gradient energy. 


100.5 
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Figure 7.4: Normalized concentration © versus normalized radial coordinate r/Ro at Cayg & 0.5 
during insertion for different particle sizes of LiFPO [96]. 


We now investigate to which extent mechanics affects the particle size depen- 
dent miscibility gap for the example of based LiFPO with E = 1 for cayg = 
(ca +cg)/2 ~ 0.5. In Fig. 7.5, we can find that the threshold value of the 
particle radius for the size independent miscibility gap expands to 43 nm in 


the mechanically coupled diffusion case, and the critical particle radius below 
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which phase segregation is inhibited increases to 9 nm. Since the strain mis- 
match between regions of different phase state induces the additional elastic 
coupling energy which scales with the width of the miscibility gap through 
the partial molar volume, minimization of the total energy leads to a reduction 
of the miscibility gap when mechanics is accounted for. The finding of a re- 
duced miscibility gap in the presence of mechanics is also consistent with the 
distribution of hydrostatic stress in Fig. 7.6 (a). There, the gradient of hydro- 
static stress induces a mechanical contribution to the driving force for diffusion 
which points to the particle center, i.e. from the high concentration phase to- 
wards the low concentration phase. As a result, the hydrostatic stress assists 
species transport towards the particle center and the concentration fields con- 
sidering the mechanics are more diffuse than that of the pure diffusion case, as 
shown in Fig. 7.6 (b). Fig. 7.6 (b) confirms that for the same particle radius 
not only the miscibility gap gets smaller, but also the width of the interface is 
widened when accounting for mechanics. This means that a comparably larger 
interface size needs to be accommodated inside the particle. This leads to a 
larger threshold value of the particle radius for the size independent miscibility 
gap and it leads to a larger particle size below which no phase segregation is 
obtained at all. On the other hand, for the small particle size the difference in 
concentration distribution between the pure and mechanically coupled diffu- 


sion cases increases as the particle size decreases. 
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Figure 7.5: Solubility limits as a function of particle 
pure and mechanically coupled diffusion 


mechanically coupled diffusion case [96]. 
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Figure 7.6: (a) Hydrostatic stress Ty versus normalized radial coordinate R/Ro at Cayg 7: 0.5 dur- 


ing insertion for different particle sizes of LiFPO with E — 1. (b) Normalized concen- 


tration versus normalized radial coordinate at Cavg © 0.5 during insertion for different 


particle sizes of LiFPO, where the solid lines represent the mechanically coupled dif- 


fusion case of E — 1 and the dashed lines 
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represent the pure diffusion case [96]. 


7.4 Evolution of the miscibility gap during insertion 


7.4 Evolution of the miscibility gap during 
insertion 


As discussed above the size dependence of the miscibility gap has been inves- 
tigated so far at a stage of average concentration of Cayg = (Ca -- cg)/2. This 
value is of particular importance for the critical particle radius below which 
phase segregation is inhibited. We now investigate how for a fixed particle size 
the miscibility gap evolves during the process of insertion. Again, we first ex- 
clude the effect of mechanics. 

Figs. 7.7 - 7.9 show that for the smaller particle sizes the solubility limits sig- 
nificantly depend on the average concentration for all three cathode materials. 
As the particle size decreases, an extended solid solution range is observed 
even in the metastable nucleation and the unstable spinodal regions of average 
concentration. In general, the miscibility gap expands as the average concen- 
tration increases. For larger particles, the system eventually enters a constant 
miscibility gap state. Interestingly, for Ro — 12.5 nm the miscibility gap of 
LMO in the whole range of phase segregated states is even smaller than that of 
NaFPO, while for larger particle sizes it is the other way around. This means 
that the suppression of phase segregation by reducing the particle size is more 
easily achieved in LMO. This is also consistent with the fact that, among the 
three materials, LMO has the largest critical radius below which phase segre- 
gation is inhibited. On the other hand, once the particle radius reaches about 
125 nm, LMO reveals a miscibility gap which is constant in the whole range 
of phase segregated states, i.e. a miscibility gap indepentdent of average con- 
centration. Compared to LMO, the threshold values of the particle radius for 
the miscibility gap to be independent of average concentration are smaller for 
LiFPO and NaFPO, namely 73 nm and 75 nm, respectively. 
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7 Particle size and average concentration dependent miscibility gap 
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Figure 7.7: (a) Solubility limits as a function of average concentration Cavg during insertion for 
different particle sizes of LMO. The red and black curves represent cq and cg, respec- 
tively. (b) Normalized concentration © versus normalized radial coordinate r/Ro for 
different time instants during insertion for three different particle sizes of LMO [96]. 
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Figure 7.8: 


(a) Solubility limits as a function of average concentration Cayg during insertion for 
different particle sizes of LiFPO. The red and black curves represent cg and cg, re- 
spectively. (b) Normalized concentration © versus normalized radial coordinate r/Ro 
for different time instants during insertion for three different particle sizes of LiFPO 


[96]. 
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7 Particle size and average concentration dependent miscibility gap 
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Figure 7.9: (a) Solubility limits as a function of average concentration Cavg during insertion for 
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different particle sizes of NaFPO. The red and black curves represent cq and cg, re- 
spectively. (b) Normalized concentration © versus normalized radial coordinate r/Ro 
for different time instants during insertion for three different particle sizes of NaFPO 


[96]. 


7.4 Evolution of the miscibility gap during insertion 


The evolution of the miscibility gap with average concentration for the smaller 
particle sizes can be explained by the minimization of system free energy for 
the example of LiFPO, as shown in Fig. 7.10. For Ro — 100 nm, the system 
reaches the maximum miscibility gap immediately once phase segregation is 
initiated (c;,4 = 0.13 for Ro = 100 nm), such that the system exhibits a misci- 
bility gap, which is independent of average concentration. For smaller particle 
sizes, for example, Ro — 12.5 nm, the system shows a dynamically increasing 
miscibility gap behavior until Cgyg = 0.49. What is the reason behind that? 
First, the interface volume is the largest at the beginning of phase segregation 
due to the shrinking core dynamics, leading to a high cost of average gradient 
energy Be, = fg 1/2 à/RÈ IVa av /V). Second, there exists a relatively 
large average gradient energy for smaller particle sizes, which is mainly con- 
trolled by the ratio A / R$. Therefore, if the maximum miscibility gap immedi- 
ately occurs once phase segregation is initiated for smaller particle sizes, the 
cost of the total system free energy would be even higher than the normalized 
chemical potential. In order to minimize the system free energy, the system 
shows a dynamically increasing miscibility gap behavior for the smaller parti- 


cle sizes before the system reaches the maximum miscibility gap. 
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7 Particle size and average concentration dependent miscibility gap 
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Figure 7.10: Normalized average system free energy Fae and normalized homogeneous free en- 
mwp 


ergy density W as function of Cayg and c, respectively, for different particle sizes 


of LiFPO during insertion. 


Now we investigate how mechanics influences the average concentration de- 
pendent miscibility gap based on LiFPO. As expected, the miscibility gap of 
the mechanically coupled diffusion case shrinks due to the suppressing effect 
of the coupling energy, as shown in Fig. 7.11. When accounting for mechan- 
ics, the initiation of phase segregation is postponed and the miscibility gap only 
exists for a smaller range of the average concentration. For the smaller parti- 
cle radius, for example, Ro — 12.5 nm, the miscibility gap of the mechanically 
coupled diffusion case always expands in almost the whole range of phase seg- 
regated states, while the miscibility gap of the pure diffusion case becomes 
constant once Cayg approaches 0.49. On the other hand, the threshold value 


for the miscibility gap to be completely independent of average concentration 
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7.4 Evolution of the miscibility gap during insertion 


expands to 100 nm in the mechanically coupled diffusion case compared to 73 


nmof the pure diffusion case shown in Fig. 7.8 a. 
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Figure 7.11: Solubility limits as a function of average concentration Cavg during insertion for dif- 
ferent particle sizes of LiFPO. The solid lines represent the mechanically coupled 
diffusion case of E — 1 and the dashed lines represent the pure diffusion case [96]. 
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8 Nonlocal species 
concentration model 


In chapter 4, a phase-field theory named NSC theory in terms of diffusion 
and phase changes has been derived. This theory incorporates two second- 
order partial differential equations involving second-order spatial derivatives 
of species concentration and an additional variable called nonlocal species 
concentration. In this chapter, we implemented the NSC theory in COMSOL 
Multiphysics® for a spherically symmetric boundary value problem of lithium 
insertion into a LMO cathode material particle of a LIB, and present the re- 
sults for spherical LMO electrode particles based on the NSC theory (see also 
Zhang and Kamlah [98]). We first introduce the spherically symmetric bound- 
ary value problem for this theory. Then, the material parameters controlling 
the interface are determined for LMO and the interface evolution is studied. 
Comparison to the Cahn-Hilliard theory shows that nonlocal species concen- 
tration theory is superior when simulating problems where the dimensions of 
the microstructure such as phase boundaries are of the same order of magni- 
tude as the problem size. This is typically the case in nanosized particles of 
phase separating electrode materials. For example, the nonlocality of nonlocal 
species concentration theory turns out to make the interface of the local con- 
centration field thinner than in Cahn-Hilliard theory. Finally, the influence of 
interface related parameters, including the penalty energy coefficient D and the 


characteristic interface length scale /, is investigated. 
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8 Nonlocal species concentration model 


8.1 Spherically symmetric boundary value 


problem 


We consider a spherical particle of radius Ro under spherically symmetric 


boundary conditions. Here, spherical coordinates are introduced and all fields 


are expressed as a function of time t and the radial coordinate 0 € r < Ro in the 


form 


(8.1) 
(8.2) 
r= Ro 
<_ 
J- n = const 
gadê i = 0 


Figure 8.1: Boundary and symmetry conditions for the spherically symmetric boundary value 


problem of the NSC theory. 


For the spherically symmetric boundary value problem, based on Equation 


(4.54), the dimensionless diffusion equation for the NSC theory can be ob- 


tained as 


—P'(1-- cl ae 
or 
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8.1 Spherically symmetric boundary value problem 


where the nonlocal species concentration is calculated by solving the dimen- 
sionless Helmholtz equation 


A a 190€ oz a 
gU 3, t7 B(c— 6) — 0. (8.4) 


The boundary conditions for the NSC theory are sketched in Fig. 8.1. At the 
surface, except for the natural boundary condition (4.47), a spatially indepen- 
dent mass flux is chosen, which has been introduced in Equation (5.30). 


In addition, at the particle center, the boundary conditions 


“0 =0, (8.5) 
on =0 (8.6) 


are imposed. Here, Equations (8.5) and (8.6) are needed to ensure the spherical 
symmetry. What is more, these two boundary conditions ensure the physical 
requirement that the flux at the particle center vanishes (see Appendix A.4). 


Finally, the initial conditions are given by 


c(r,0) =0, (8.7) 
é(r,0) =0. (8.8) 


The cathode material LMO is chosen for the simulations. We consider a typical 
value of Ro = lum for the particle radius. The material parameters for the 
cathode material LMO have been summarized in Table 7.1, and the parameters 


B and / will be discussed in the following section. 
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8 Nonlocal species concentration model 


8.2 Interface evolution 


Here, we will consider quasistatic insertion of lithium into a particle of LMO 
at C — 0.001. First, we consider six different cases of p and / with the same 
value of the product A = BI? = 7 x 10-15 m?. Figs. 8.2 and 8.3 show the nor- 
malized lithium concentration © and the normalized nonlocal concentration © 
along the radial coordinate at a state of average concentration of Cayg = 0.5 for 
different cases of p and /. In Fig. 8.2 , we find that the lithium concentration 
plots of different cases are almost the same except for the first case (B = 0.001 
and / — 83.7 nm) which shows a homogeneous state. All other cases exhibit 
an extremely small difference in the interface situated between the two phases 
as shown in the inset. Also, Fig. 8.3 shows almost the same nonlocal con- 
centration plots for the different cases. As illustrated in Fig. 8.4, for the last 
case (B = 1000 and / = 0.084 nm), the lithium concentration is the same as the 
nonlocal concentration, both of which are also identical to the lithium concen- 
tration from the Cahn-Hilliard model. According to Eq. (5.46), for the limit 
case l — 0 the nonlocal concentration is equal to the lithium concentration c, 
meaning that the Cahn-Hilliard model is obtained and the full nonlocal effect 
is reduced to a weak one. Therefore, p — 1000 in our simulations is enough to 
approach the results from the Cahn-Hilliard theory. This is also confirmed by 
Di Leo et al. [94]. 
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Figure 8.2: Normalized lithium concentration © versus normalized radial coordinate r/Ro at Cayg = 


0.5 during insertion for different cases of B and / [98]. 
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Figure 8.3: Normalized nonlocal concentration © versus normalized radial coordinate r/Ry at 


Cavg = 0.5 during insertion for different cases of B and / [98]. 
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Figure 8.4: Comparison of the concentration plots at Cayg = 0.5 during insertion from the NSC 
theory and the Cahn-Hilliard theory, where B = 1000, / = 0.084 nm, and A — 7 x 10718 
2 
m^. [98] 


Now we focus on the interface evolution shown in the insets of Figs. 8.2 and 
8.3. As I increases, the interface of the lithium concentration field becomes 
steeper, but the nonlocal concentration field gets more diffuse, which means 
that the spatial interaction volume which contributes to the nonlocal effect is 
larger. The reason is that due to the fixed product A, the increase of the spatial 
interaction volume is compensated by a shrinking of the interface which is 
governed by the characteristic interface length scale / [113]. Therefore, the 
nonlocality makes the interface of the local concentration field c steeper. In 
addition, we should notice that the first case (B = 0.001 and / = 83.7 nm) 
shows no phase segregation, which is explained by Fig. 8.5. Here, zero system 
penalty energy corresponds to a homogeneous system. For the fixed product 
A, as | increases, the penalty energy increases. When / reaches 83.7 nm, the 
cost of the penalty energy \PP"alty at a phase segregated state would be too 


high so that as a result of the competition between the different contributions 
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8.2 Interface evolution 


to the total system energy the system stays in the homogenous state in order to 
minimize the total system free energy. Therefore, the penalty energy vanishes 
in this case, and it can be seen that the contribution of the penalty energy to the 


total system free energy can lead to suppression of phase segregation. 
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Figure 8.5: Normalized penalty energy ‘¥?¢"“""" versus average concentration Cavg during insertion 
for different cases of B and / [98]. 


Next, we discuss the evolution of the interface thickness d of the lithium con- 
centration field, the interface thickness dnon of the nonlocal concentration field, 
and the characteristic interface length scale / with increasing D during inser- 
tion, as shown in Fig. 8.6. Here, d or dnon is defined as the distance over 
which the concentration deviates 0.05 from the two binodal concentrations of 
the lithium concentration field or the nonlocal concentration field [113, 141]. 
We find that d increases with increasing B, while both dnon and / decrease. As 


p increases, d is more close to dnon, and they are the same when ß reaches 
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1000 as shown in the inset, when approaching the solution from the Cahn- 
Hilliard theory as discussed before. However, the difference |d — /| reduces as 
B decreases. When ß is 1.5, / and d are equal to 2.16 nm and 5.66 nm, re- 
spectively, and thus / is on the order of magnitude of d. Therefore, in view of 
the interpretation of / as the characteristic interface length scale, we can regard 
the two values 1.5 and 2.16 nm as reasonable rough estimates of B and / for 
the cathode material LMO, respectively, when fixing A = 7 x 107!8 m?, which 
leads to a lower interface thickness of 5.66 nm in the NSC theory compared to 


the interface thickness of 8 nm in the Cahn-Hilliard theory. 


NSC, C=0.001, c... -0.5, R,=1 um 


497x107? m? 


Interface thickness parameters (nm) 


Figure 8.6: Evolution of the interface thickness d of the lithium concentration field, the interface 
thickness dnon of the nonlocal concentration field, and the characteristic interface length 
scale / at cayg = 0.5 according to the NSC theory with A = 7 x 10-73 m? during inser- 
tion [98]. 
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8.3 Comparison with Cahn-Hilliard model 


8.3 Comparison with Cahn-Hilliard model 


As the truly nonlocal NSC theory is an extension of the Cahn-Hilliard theory, 
we will now compare the results from the two theories. First, we consider the 
electrode particle as being microsized, say, Ro = 1 um. Then, as illustrated 
in Fig. 8.7, the lithium concentration plots obtained from the two theories 
are almost the same. This means that, for a microsized electrode particle, the 
nonlocal effect is very weak, and the Cahn-Hilliard theory yields the almost 
same results as the NSC theory. Therefore, the Cahn-Hilliard theory is well 
capable of describing the motion of the sharp interface between the two phases 
for a microsized electrode particle. However, when the electrode particle radius 
is considered as being nanosized, say, Ro = 100 nm, we find in Fig. 8.8 that the 
lithium concentration at the interface between the two phases obtained from the 
Cahn-Hilliard theory is different from that of the NSC theory. The interface of 
the lithium concentration field c considering the truly nonlocal effect is steeper 
than that from the Cahn-Hilliard theory, and the nonlocal concentration field 
€ is the most diffuse among the three concentration plots. Indeed, the ratio 
of the characteristic interface length scale to the particle radius //Ro in the 
nanoparticle is larger than that in the microparticle, and, as a result, the larger 
ratio of the spatial interaction volume to the particle volume contributes to a 


stronger nonlocal effect. 
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Figure 8.7: Comparison of the concentration plots at Cayg = 0.5 during insertion from the NSC 
theory and the Cahn-Hilliard theory, where Ro = 1 um, p = 1.5, 1 — 2.16 nm, and 
A =7 x 10738 m? [98]. 
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Figure 8.8: Comparison of the concentration plots at cayg = 0.5 during insertion from the NSC 
theory and the Cahn-Hilliard theory, where Ro = 100 nm, B = 1.5, 1 = 2.16 nm, and 
A=7x 10738 m? [98]. 
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8.3 Comparison with Cahn-Hilliard model 


Next, we consider a microparticle with a large microstructure scale. To this 
end, we choose A = 7 x 107? m?, to vary for the sake of comparison the 
interface thickness in both, NSC and Cahn-Hilliard theory. For fixed p = 1:5; 
we obtain / = 68.31 nm. Fig. 8.9 shows the results from the two theories. It 
is found that the lithium concentration from the Cahn-Hilliard theory is clearly 
different from that of the NSC theory. The lithium concentration field from the 
Cahn-Hilliard theory is more diffuse than that from the NSC theory, and the 
interface of the nonlocal concentration field is very thick, which means that the 
spatial interaction volume is very large. This can be explained by the increase 
of the characteristic interface length scale / enlarging the spatial interaction 


volume. As a result, the nonlocal effect is very strong. 
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Figure 8.9: Comparison of the concentration plots at cayg = 0.5 during insertion from the NSC 


theory and the Cahn-Hilliard theory, where Ro — 1 um, B — 1.5, 1 — 68.31 nm, and 
A =7 x 107" n? [98]. 
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8 Nonlocal species concentration model 


8.4 Influence of the penalty energy 
coefficient p 


As the interface thickness d is determined by the two parameters B and / in 
the NSC theory, we will study their influence. First, we vary the value of 
B but fix 1 (| = 2.16 nm). Figs. 8.10 and 8.11 show the effect of B on the 
normalized lithium concentration c and the normalized nonlocal concentration 
€ at a state of average concentration of Cayg = 0.5. We find that both, d and dnon, 
increase as B increases, but a homogeneous state is obtained for the last case 
(B — 15000). This can be explained with the help of Figs. 8.12 and 8.13. Fig. 
8.12 shows that with increasing values of ß, the interface thicknesses d and 
dnon of © and €, respectively, tend towards each other asymptotically. This in 
turn leads to the plots of c and é to coincide more and more and, by the way, to 
approaching the solution of the Cahn-Hilliard theory. The results show that the 
term (c — &/ dominates the trend of increasing ß in the penalty energy density 
(4.51). As a result, the system penalty energy in Fig. 8.13 is gradually reduced 
with increasing value of B. Finally, when f reaches 15000, the system penalty 
energy vanishes, and a homogeneous state occurs. As discussed before, the 
asymptotic coincidence of the NSC theory with the Cahn-Hilliard theory for 


large values of ß means a likewise asymptotic reduction to weak nonlocality. 
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Figure 8.10: Normalized lithium concentration © versus normalized radial coordinate r/Ro at 


Cayg = 0.5 during insertion for different B, where l = 2.16 nm [98]. 
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Figure 8.11: Normalized nonlocal concentration © versus normalized radial coordinate r/Ro at 


Cayg = 0.5 during insertion for different B. where / = 2.16 nm [98]. 


135 


8 Nonlocal species concentration model 


2 21 p 214 2143 


M 
a 
o 


200 


a 
o 


o 
o 


Interface thickness parameters (nm) 
o 
3 


1E-14 


Figure 8.12: Evolution of the interface thickness d of the lithium concentration field, the inter- 
face thickness dnon of the nonlocal concentration field and the characteristic interface 
length scale / at Cavg = 0.5 during insertion, where / = 2.16 nm [98]. 
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Figure 8.13: Normalized penalty energy PPenalty versus average concentration Cavg during inser- 


tion for different B. where / = 2.16 nm [98]. 
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8.5 Influence of the characteristic interface 
length scale / 


Next, we vary the value of / but fix B = 1.5. Figs. 8.14 and 8.15 show the effect 
of 7 on the normalized lithium concentration c and the normalized nonlocal 
concentration ĉ at a state of average concentration of Cayg = 0.5. As it has to 
be expected, increasing the characteristic interface length scale / broadens d 
and dnon, and increases the difference between é and 6, which indicates that 
the nonlocal effect is enhanced. This is consistent with the increase of the 
maxima of the system penalty energy in Fig. 8.16 as / increases. Interestingly, 
comparing the case (B = 1.5, A = 7 x 10 ^ m?) in Fig. 8.14 to the case (l = 
2.16 nm, A = 7 x 107^ m?) in Fig. 8.10, we find that it is a phase segregated 
state for the former case but a homogeneous state for the latter case, although 
the product A is the same. Indeed, ß in the latter case is far more than 1000, 
which means that the result from the Cahn-Hilliard theory is approached. On 
the other hand, the former case represents a result reasonably accounting for 
a truly nonlocal effect due to the fact that the characteristic interface length 
scale / is on the order of magnitude of the interface thickness d. Therefore 
we conclude that, the Cahn-Hilliard theory is not capable of describing the 
interface evolution of a microstructure possessing a length scale comparable 
to the total system size. Now, we focus on the evolution of d and dnon with 
increasing / for fixed D, as shown in Fig. 8.17. It is found that the difference 
|d — dnon| initially increases as / increases. Beyond / = 68.3 nm, the difference 
|d — dnon| gets smaller. Finally, for the case of / = 683 nm and A = 7 x 107! 
m2, which is not shown in the figures, there is a homogeneous state and, as a 
consequence both d and dnon vanish with the trivial consequence |d — dnon| = 0. 
The reason for this change of the trend of |d — dyon| lies in the competition 
of the different terms contributing to the system free energy. As the system 
penalty energy is about to get too large (see Fig. 8.16), a homogenous solution 


yields the minimum system free energy. In other words, too large a penalty 
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energy suppresses phase segregation. In addition, d is always on the order of 
magnitude of / in the different cases. Therefore, the value 1.5 can be confirmed 


again as a reasonable rough estimate of B for the NSC theory. 
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Figure 8.14: Normalized lithium concentration © versus normalized radial coordinate r/Ro at 
Cave = 0.5 during insertion for different /, where B = 1.5 [98]. 
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Figure 8.15: Normalized nonlocal concentration ĉ versus normalized radial coordinate r/Ro at 
Cave = 0.5 during insertion for different /, where B = 1.5 [98]. 
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Figure 8.16: Normalized penalty energy \PPeralty versus average concentration Cave during inser- 
tion for different /, where B = 1.5 [98]. 
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Figure 8.17: Evolution of the interface thickness d of the lithium concentration field, the interface 
thickness dnon of the nonlocal concentration field and the interfacial characteristic 
interface length scale / at Cavg = 0.5 during insertion, where B = 1.5 [98]. 
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The occurence of phase segregation within electrode particles, which is con- 
sidered as a resource of large mechanical stresses, can be a possible origin for 
the battery capacity fade. Thus, it is the objective of this work to study phase 
changes and mechanical stresses in electrode particles by means of phase- 
field simulations. This work was first dedicated to the understanding of phase 
changes, mechanical deformation, and stress evolution in NaFPO electrode 
particles of NIBs. Furthermore, the particle size and average concentration 
dependent miscibility gap for the three cathode materials LMO, LiFPO, and 
NaPPO were studied. Finally, a NSC theory, which is an extension of the 
Cahn-Hilliard theory, was introduced. 

We formulated a non-linear Cahn-Hilliard type model with mechanics. For 
the mechanical part, besides SDT two different finite deformation elasticity 
formulations were introduced and compared. This work has addressed the 
theoretical formulation, numerical implementation, and application of the me- 
chanically coupled phase-field theory. 

To start with, a phase-field model for the cathode material NaFPO of NIBs 
was studied for the first time, using a phase-field theory coupling the Cahn- 
Hilliard equation to finite deformation elasticity. As a major novelty, the ma- 
terial parameters for NaFPO were determined. For example, 0, 05, and A, 
all of which are the key parameters in the phase-field model, were determined. 
The determination of these key parameters provides a significant input for the 
future phase-field work for NaFPO. We implemented the fourth-order nonlin- 
ear initial-boundary-value problem of the model in COMSOL Multiphysics® 


to solve by the finite element method the spherically symmetric problem of 
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sodium insertion into or extraction from a NaFPO particle of NIBs. Our model 
captures the important feature that distinguishes NaFPO from LiFPO, i.e, 
phase segregation into a sodium-poor phase FePO, and a sodium-rich phase 
Nay/3FePO4. The difference for the concentration and stress plots between 
SDT and the finite deformation theories can not be neglected for NaFPO. In 
particular, the stresses in the two phases are higher for SDT. The above differ- 
ence suggests that the finite deformation elasticity is preferred for the cathode 
material NaFPO. This is different from that SDT has a sufficient capacity to 
represent the deformation of LiFPO. Although the results from the two finite 
deformation theories are almost the same, according to our experience, elas- 
ticity based on logarithmic elastic strain is numerically more efficient than 
elasticity based on elastic Green strain. On the other hand, we found that, 
for the same stage of average concentration, the miscibility gap of NaFPO is 
smaller than that of LiFPO, but the stresses in a NaFPO particle during phase 
segregation are larger, which may explain the existence of a wide range of 
solid solution Na,FePO for 2/3 « x « 1 to avoid even higher stresses in the 
NaFPO particle. In addition, the suppression of phase segregation by the elas- 
tic strain energy is more pronounced in NaFPO compared to LiFPO. Here, we 
just focus on phase changes in the two-phase region of NaFPO (0 « x « 2/3), 
ignoring the single-phase region (2/3 « x « 1). As an outlook, the current 
model may be extended into the whole region (0 « x « 1). As both the partial 
molar volume and Young's modulus are regarded as constants here, another 
major subject of future study is the effect of the concentration dependence of 
these two quantities for NaFPO. 

For the second aim of this work, by means of a phase-field theory coupling the 
Cahn-Hilliard equation to finite deformation elasticity, we studied the depen- 
dence of the miscibility gap on particle size and average concentration for the 
three cathode materials LMO, LiFPO, and NaFPO. We found that, a reduced 
miscibility gap in the center of the spinodal region of average concentration 
is found for radii below 55.5 nm, 20.5 nm, and 17.5 nm for the three cathode 
materials LMO, LiFPO, and NaFPO, respectively. The miscibility gap shrinks 
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to accommodate in particles of reduced size the gradient of concentration in 
the interface between phases. The critical radius below which phase segrega- 
tion is inhibited for LMO is R. = 10.9 nm, which is larger than the critical 
radii of 6 nm for LiFPO and 6.4 nm for NaFPO. Concerning the evolution of 
the miscibility gap during the insertion process, it was found that for smaller 
particles it increases as average concentration increases. However, the mis- 
cibility gap is constant during the whole process of insertion for radii above 
125 nm, 73 nm, and 75 nm for the three cathode materials LMO, LiFPO, and 
NaFPO, respectively. The average concentration dependent miscibility gap is 
physically explained by the minimization of system free energy. On the other 
hand, when mechanics is accounted for, minimization of the total energy leads 
to a reduction of the miscibility gap and a widened interface region. Among 
the three investigated cathode materials, mechanics has the strongest influence 
on LMO on the miscibility gap, while the tensile stresses in this material are 
the smallest. The elastic strain energy leads to a larger threshold value of the 
particle radius for the size independent miscibility gap, and a larger critical 
particle size below which phase segregation is completely inhibited. Our work 
suggests that one way to improve the mechanical stability of cathode mate- 
rials is that the miscibility gap is suppressed by choosing sufficiently small 
particle sizes according to the respective material. For the future challenge 
to design nanoscale insertion materials, LMO possesses the largest particle 
size for expanding the solid solution range in which particle stresses are elim- 
inated among the three investigated cathode materials. Here, neither surface 
tension nor surface wetting [50, 51, 142] is taken into account. According to 
Stein et al. [142], surface tension plays an increasingly unnegligible role in the 
electro-chemo-mechanical behavior of electrode nanoparticles as the particle 
size decreases. Investigating the effect of surface tension on the miscibility 
gap in phase separating insertion materials can be subject of future work. In 
order to avoid computationally expensive three-dimensional simulations, the 
cathodic particles are assumed to be of spherical symmetry under spherically 


symmetric boundary conditions. However, the particle shape has an important 
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influence on the phase segregation and the stress state [88, 95, 116]. According 
to Zhao et al. [95], the core-shell phase segregation may not occur due to the 
geometric anisotropy in plate- and needle-like particles, and the stress magni- 
tudes in the particles of the above two types can be different compared to those 
in particles of spherical symmetry. Investigating the effect of the particle shape 
may also be subject of future work. 

Finally, a NSC theory was introduced from a general formulation of a non- 
local free energy density by using the Green's function as weight function to 
describe diffusion and phase changes in a material showing phase segregation. 
This theory was shown to be a generalization of the classical Cahn-Hilliard 
theory. In contrast to the Cahn-Hilliard theory which is weakly nonlocal, the 
NSC theory is truly nonlocal and is based not only on the species concentra- 
tion c but also on the nonlocal species concentration ¢. The NSC theory is 
governed by two coupled second-order PDEs which are amenable to being 
solved with a standard finite element implementation and, thus, is computa- 
tionally less demanding than the fourth-order Cahn-Hilliard equation. In the 
NSC theory, the nonlocal free energy density is split into two parts, which are 
the penalty energy density y^*"*!''* and the variance energy density w’iance, 
The interface thickness d is controlled by two independent parameters, namely 
the penalty energy coefficient p and the characteristic interface length scale /. 
We implemented the NSC theory in COMSOL Multiphysics® to solve by the 
finite element method the spherically symmetric problem of lithium insertion 
into a LMO particle of a lithium ion battery. It was found that ß = 1.5 and 
l — 2.16 nm can be regarded as reasonable rough estimates for the cathode ma- 
terial LMO. The nonlocality makes the interface of the local concentration field 
steeper, when compared to the Cahn-Hilliard theory. Also, the NSC theory is 
more accurate than the Cahn-Hilliard theory for a nanoparticle of an electrode 
material or, more general, a comparable scale of the microstructure in relation 
to the system size. A high system penalty energy can lead to suppression of 


phase segregation. On the other hand, both interface thicknesses d and dnon 
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increase as D or l increase. Furthermore, increasing B weakens the nonlocal 


effect while increasing / enhances the nonlocality. 
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A Appendix 


A.1 Helmholtz equation (4.46) as governing 
partial differential equation for 
nonlocal species concentration 


The Green's function is the solution of the Helmholtz equation (4.46) with the 


source term replaced by the Dirac function Ó (x — y), i.e., the solution of 
G:F) -P VGE) = 8(X— y), (A.1.1) 


for the boundary condition (4.47). Integrating Equation (A.1.1) over the prob- 


lem domain ZZ with respect to spatial position y, one obtains 


/ G(yV = / 8( Nav +22 f V?G(;5)dV 
B B B 


2 1+P | VG(:y) -HAA 
OB 
= A (A.1.2) 


where the term related to /* vanishes by the divergence theorem and the bound- 
ary condition (4.47). Consequently, if the Greens’s function G(X; y) is contin- 
ued by zero values outside the problem domain Z, it satisfies the natural weight 
function property (4.28), and, thus it is an admissible choice of the weight 
function, ie. @();X) = G(J;X). Notice that G(¥;x) = G(X;y) since Green's 


functions are symmetric for linear problems with constant coefficients [114]. 
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Thus, we introduce the nonlocal species concentration according to Equation 
(4.35) as 


I 


A 
(x) = eT Lnso:nav. (A.1.3) 


In view of the choice of G as the weight function, the definition (4.36) together 
with Equation (A.1.2) yields 


A A 
p- n = —, (A.1.4) 


Cmax 


such that Equation (A.1.3) simplifies to 


éx) = J E005. (A.1.5) 


By the definition of the Dirac function 


B ss L SEa (A.1.6) 
holds. Now, substituting Equation (A.1.1) into Equation (A.1.6) gives 
eg) = [ EOE- | PEVE) 
e)G:3)av — P? | EC 


B 
&(x) - PV?8(X). (A.1.7) 


Il 
T 


Therefore, the inhomogenous Helmholtz equation (4.46) is the governing equa- 
tion for nonlocal species concentration when is defined according to the in- 
tegral representation (4.35) with the Green's function as the weight function. 


As an additional implication, as also shown in [114], the coefficients related to 
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A.2 Derivation of the system free energy in the nonlocal species ... 


derivatives of fourth-order and higher of c in Equation (4.45) need to vanish, 
such that equation (4.45) simplifies to the Helmholtz equation (4.46). 


A.2 Derivation of the system free energy in 
the nonlocal species concentration 
theory 


The free energy density (4.27) in the nonlocal model is reformulated as 


2 


wG) = VPE) +54 f oP EE-E) av 


1 </> dre (— A AMA 
+54 | ot) - 26)" — ep) (2) - 8D 'av. 
(A.2.1) 
Note that y is the variable of spatial integration. Using Equation (4.36), the 
second term on the right hand side of Equation (A.2.1) gives 


1 


A f OEE- AV = zenb (EE-E) 


= wrenalty (E(x) : é(X) ) : (A.2.2) 


Using Equations (4.35) and (A.1.2), the third term on the right hand side of 
Equation (A.2.1) then yields 


La f PEE- -EE-E av 


= (A — Cmax EEEE) — -AG(X)? + A s e(p)c(y)dV. | (A2.3) 
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The first term on the right hand side of Equation (A.2.3) vanishes by Equation 
(A.1.4), and Equation (A.2.3) then can be manipulated further: 


- 3^8 +54 Ph w(p)é(3t)2aV 
i ys oX(p)e()^av — &()?) 
= T TES o(p)e(y)?dv + &(x)* — 26(x)?) 
= T (d, a(p)e(y)^av 4-é(x?? — 2&(x) [e (eX av) 
z an olp) (E0) -Eav 
Lo yrariance DN (A.2.4) 


As shown in appendix A.1, the Green's function G(y;X) of the Helmholtz equa- 
tion (4.46) satisfies the natural weight function property (4.28). Thus, it is 
chosen as the weight function of the NSC theory. As a result, we obtain the 
system free energy (4.50) in which the nonlocal free energy density is com- 
posed of two terms, namely the penalty energy density yP“*@l'y and variance 


energy rannte : 
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A.3 Derivation of the final coupled diffusion equation ... 


A.3 Derivation of the final coupled 
diffusion equation of the spherically 
symmetric boundary value problem 


A.3.1 The coupled Cahn-Hilliard equation with 
elastic Green strain 


Since Grad iz" (cg, E^) is extremely complicated and long, it makes sense to 


use the following substitutions: 


s = Q(cg — co), (A.3.1) 
q = x (A.3.2) 
w = E (A.3.3) 
V 
n = Du (A.3.4) 
X = 0495, (A.3.5) 
Y = (49), (A.3.6) 
Z = (1+9, (A.3.7) 
Q = (149), (A.3.8) 
W = (1+w). (A.3.9) 


According to Table 5.1, all the above introduced parameters (A.3.1 - A.3.9) are 
dimensionless. 
For the gradient of the chemical potential, Grad 17"? and Grad a have been 


expressed in Equations (5.25) and (5.26), respectively. Again, here c and r in 
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Equations (5.25) and (5.26) should be respectively replaced by cr and R. Using 
Equation (5.80), Grad L7? is expressed as 


20 724 
3 SZ OY 20-92 IR 


Graduf'(cg,E^) = GQ „(20 1)( ) 


—)(1+w))(WZ-1)) 


+2Z(1+w)(4— 7.) | ë 


42n(- or = wy) 
{—1 A ER 


1 (A.3.10) 


A.3 Derivation of the final coupled diffusion equation ... 


Substituting Equations (5.25), (5.26), and (A.3.10) into Equation (5.68), we 


can derive the final dimensionless coupled diffusion equation 


-2 0€ 0 E C 
ZSSR | Hees = Z 
R JF t IR | R (1 t OoCg(1 €g)) = +R 
Por à 2 OPER 2 ACR 
ƏR? RƏR? R2 OR 


a= [1 2 Os dq 2 Os 
[cà 520 1)( 3 aR et t2(l+a)Z a5) + ( 3IW3g 


M2z($ - )(1+w))(WZ-1)+2n(-1+ se +QZ)+WZ) 


a a 
Z EZü og) — £wy £5. e 2Z( 4-wy( 


)) *2n( 


Nee WC eon ut ) 


3 3 
Os | dq 2 Os q w 
OY = Ara) RZ aag t2wt N) — wz] 


ox 55 (l+q)=sY+ 


Ep: 


(ZQ —1)QY WZ-—1)WWY 


QY uni 1+5( 1+02)+wz)| |. 


(A3.11) 
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A.3.2 The coupled Cahn-Hilliard equation with 
logarithmic elastic strain 


Since Grad u£" (cr, E 


fog) is extremely complicated and long as well, it makes 


sense to use the following substitutions: 


x* = pay (A.3.12) 
Y* = (+83, (A.3.13) 
7 2, (A.3.14) 
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A.3 Derivation of the final coupled diffusion equation ... 


According to Table 5.1, all the above introduced parameters (A.3.12 - A.3.14) 
are dimensionless. 


Using Equation (5.92), Grad JP" (cr,E},,) is expressed as 


2InÀ* / 19 a 

Grads? (or Bf) = 60/22 (Cicerone zii) 
1 

4In A? 1 ðs q w 
NE 2 __ a" y* *(12 7 7 

A: ( gon! EM, D) 
-2n (InA* +21nA?) 

1 T- E „gg 
(a 3A apt vera ) 


2 lds., «q4 
A: (3507 (1+q)+Z*( 


„9q In A7 | * 19s * i 


„og | nA’ „gg 4lnA? oq 
+Z ) "je oR ae (1+g)X ) 


Ind; ( Y'9s «4 
iiem (- Ie - 3) 
In A7 q Wx 4lnA7 , 05 

^n - c EO a) 
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Ae \ 30R 
2 los., 4,0 w 
Q . 22 * 
(cacao Eu 


+2n (InA° *-21nA7) (autor 


los * „9q 

(3337 (l+q)+Z 3) 
Q 94 y. | 4Q 
3A¢ OR | OAS 


20. 


4Q InA? 
2 (1+w)Y* +2nQ (Ind +21nA°) 


ll y*(14.9) S DIK 
B rn ia 
(A.3.15) 
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A.3 Derivation of the final coupled diffusion equation ... 


Substituting Equations (5.25), (5.26), and (A.3.15) into Equation (5.68), we 


can derive the final dimensionless coupled diffusion equation 


-, 0G a > OCR aA 
20CR PE = R 2^. A 
0 = R ET +53 | R (1 j O2€r(1 €g)) JR +R go R) 
or 2 8ER 2 dcr 2. 
(ae eor) fto 


PR 21nA* los., ,9q 
[en 5 (oes) 


4lnA? 19s " i xd w 
ur" ( aon CINP G-p) 


1 1 9 ð 
+2n (InA* +21nA?) (Er SSY +q) +25) 
1 


1 


ET (er 


, 1 E 1s, «9q 
+G(1 +s) |— 3 (is IR er (l+q)+Z a) 
Ina’ * 1 ds * , 0q 
= m ( yy er (14-3) -Z 5) 


pA ðq 4lnA* 


IE ae TAA sa) 


40 (1 wh 1ds "x i 


In A* 19 
Sa (ar ten z 6-3) 


Ae? 30R 
Ind? q w 41In A? Os 
2 i té Y* 2 1 X* d 
ie R R wm T) 
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A.4 


+2n (InA” +21nAf) | 


310° OR 
E ae ar 
aeter ( rd Ew) zw 2) 
Zi e ew) ape 
2nQ (inà? +21nA?) 
( 3,0 Hq) ||| (A.3.16) 


Lithium flux at the particle center for 
the spherically symmetric boundary 
value problem of the nonlocal species 
concentration theory 


Here, we will demonstrate that the two boundary conditions (8.5) and (8.6) en- 


sure that the flux at the particle center vanishes. This is a physical requirement 


in order to guarantee the conservation of lithium matter within the particle. 
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A.4 Lithium flux at the particle center ... 


Using Equations (4.4), (4.37), and (4.10), we can obtain the normalized lithium 


flux in the one-dimensional particle model as 


"S RT. 
ar PSN S 


Sıl Ox 


J=-(1+ē(1—2)) » (A.4.1) 


where r — r/Ro. Then, according to Equations (8.5) and (8.6), the normalized 


lithium flux at the particle center is 
F(0,t) — 0. (A.4.2) 


Therefore, the two boundary conditions (8.5) and (8.6) ensure the physical 


requirement that the lithium flux at the particle center vanishes. 
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Sodium-ion batteries become an attractive alternative to lithium-ion batteries. 
Most storage materials exhibit phase changes during operation. The respective 
phases of a storage material possess different lattice constants giving rise to a 
strain mismatch which, in turn, causes mechanical stresses and, thus, leads to 
damage of the electrode particles. 

In this work, for non-linear Cahn-Hilliard type models describing diffusion, a 
thermodynamical framework for the coupling with mechanics is presented. First, 
a phase-field model for the cathode material Na,FePO, is studied for the first time 
to understand phase changes, mechanical deformation, and stress evolution in 
Na,FePO,. Furthermore, we study the particle size and average concentration 
dependent miscibility gap of the nanoscale insertion materials Li,Mn,O,, Li,FePO,, 
and Na,FePO,. Finally, we introduce the nonlocal species concentration theory 
from a nonlocal free energy density, compare this theory to the Cahn-Hilliard 
theory, and show how the nonlocality influences the results. 
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